[R] Matrix package,solve() errors and crashes Please help

[Martin Maechler]

>>>>> "SS" == Surendar Swaminathan <surendar.swaminathan at gmail.com>
>>>>>     on Fri, 15 May 2009 15:55:23 -0700 writes:

    >> Hello All,
    >> 
    SS> Please help me with this problem.I have been having this problem for over a
    SS> month now and I could not find any information.I later realised that error
    SS> is with MATRIX package.

[...]


    SS> This is my graph object.
    SS> The file 'Bonacich Power.RData' (1.5 MB) is available for download at
    SS> <
    SS> http://dropbox.unl.edu/uploads/20090522/7a0d0313f21fd6a8/Bonacich%20Power.RData

[...]

and you have asked almost the idetnical question 4 weeks ago,
and got a first answer, but then not a *real* answer.
I'm sorry for that, since I had started answering you, 
and even more as I *did* solve the problem back then.. (April 22).

more below ...

    SS> computer Configuration

    SS> *WIndows XP service Pack 3 .0 GB RAM*

    SS> I am using SPARSE matrix to solve the problem
    SS> This is the code I use to obtain bonpower using Sparse Matrix &
    SS> alternatively the code is in the following website

    SS> http://igraph.wikidot.com/r-recipes#toc6

I'm really appaled that "R-recipes" on such a website give 
such an unprofessional advice


{the following is so grabled because of your (implicit)
 insistence on using HTML-ified e-mail .. }

    SS> *bonpow.sparse <- function(graph, nodes=V(graph), loops=FALSE,
    SS> exponent=1, rescale=TRUE, tol=1e-07) {*
    SS> *  ## remove loops if requested
    SS> *
    SS> *  ## sparse adjacency matrix
    SS> d <- get.adjacency(graph, sparse=TRUE)*
    SS> *  ## sparse identity matrix
    SS> id <- spMatrix(vcount(graph), vcount(graph),
    SS> i=1:vcount(graph), j=1:vcount(graph),
    SS> x=rep(1, vcount(graph)))
    SS> id <- as(id, "dgCMatrix")*
    SS> *  ## solve it
>>>>  SS> ev <- solve(id - exponent * d, tol=tol) %*% degree(graph, mode="all")*
    SS> *  if (rescale) {
    SS> ev <- ev/sum(ev)
    SS> } else {
    SS> ev <- ev * sqrt(vcount(graph)/sum((ev)^2))
    SS> }*
    SS> *  ev[as.numeric(nodes) + 1]
    SS> }*
    SS> **
    SS> I realised that the error is in Matrix Package
    SS> *Error: cannot allocate vector of size 3.3 Gb
    SS> In addition: Warning messages:
    SS> 1: In solve(id - exponent * d, tol = tol) :
    SS> Reached total allocation of 1535Mb: see help(memory.size)
    SS> 2: In solve(id - exponent * d, tol = tol) :
    SS> Reached total allocation of 1535Mb: see help(memory.size)
    SS> 3: In solve(id - exponent * d, tol = tol) :
    SS> Reached total allocation of 1535Mb: see help(memory.size)
    SS> 4: In solve(id - exponent * d, tol = tol) :
    SS> Reached total allocation of 1535Mb: see help(memory.size)*
    SS> **
    SS> sessionInfo()
    SS> R version 2.9.0 (2009-04-17)

     [ ............. ] {irrelevant here}


    SS> Please help Matrix Experts

Well, I've marked the one line above

>>>> ev <- solve(id - exponent * d, tol=tol) %*% degree(graph, mode="all")*

<begin abusive language; excuse in advance ..> 
which really contains stupid advice.
Yes, this is strongly put, but I think people who tell you to
solve

    A x = b  (A [n x n] matrix)

in R by	     x <- solve(A) %*% b
should not give advice on programming at all.
<end abusive language>

One of the proper advices would use
       	     x <- solve(A, b)

and indeed that's even more crucial in the case where  A  is
sparse matrix:  solve(A) is never sparse for sparse A (*),
but   solve(A, b)  can happen with a fast (and memory-efficient)
algorithm, and that's what happens also when you use the Matrix package.
I append a version of your function that is also slightly
improved in other places.

(*) and that's why your memory blows up when A is large and sparse !

    SS> Thanks in advance
    SS> Nathan

    SS> [[alternative HTML version deleted]]

(still....   It would *really* help if you used a more sensible e-mail
	     configuration  }

Regards,
Martin Maechler

-----------

Here's the corrected bonpow.sparse() function , actually with
extra code for diagnostic output, etc.
You should probably simplify it again :

## I got help from IGRAPH community to use sparse Matrix
##
## http://igraph.wikidot.com/r-recipes#toc6
##
### "Enhanced" by Martin Maechler:
bonpow.sparse <- function(graph, nodes = V(graph), loops = FALSE,
                          exponent = 1, adj.type = "both",
                          trace = TRUE,
                          rescale=FALSE, tol=1e-07)
{
    stopifnot(require("igraph"),
	      require("Matrix"))

    if(trace) {
	c.width <- 30
	C1 <- function(s) cat(sprintf("%-*s .. ", c.width, s))
	C2 <- function() cat("[Ok]\n")
    } else { C1 <- C2 <- function(...) {} }

    ## remove loops if requested
    if (!loops) {
	C1("simplify()ing graph")
	graph <- simplify(graph, remove.multiple=FALSE, remove.loops=TRUE)
	C2()
    }

    ## sparse adjacency matrix
    C1("d <- get.adjacency(., sparse)")
    d <- get.adjacency(graph, type = adj.type, sparse=TRUE); C2()

    if(trace >= 2)
	cat("class(d): ", class(d),"\n")
    if(!is.directed(graph)) {
	## MM: unfortunately  "igraph" does not return a *symmetric*
	## --- sparse matrix directly, saving space and time,
	## so we at least do it now :
	C1("d <- as(d, \"symmetricMatrix\")")
	d <- as(d, "symmetricMatrix"); C2()
    }

    ## sparse identity matrix
    vg <- vcount(graph)
    if(FALSE) { ## "non-sense" :
	C1("spMatrix(.) for Diagonal(vg)")
	id <- spMatrix(vg, vg, i=1:vg, j=1:vg, x = rep(1, vg)); C2()
	C1("	--> as(., \"dgCMatrix\")")
	id <- as(id, "dgCMatrix"); C2()
    }
    else
	id <- Diagonal(vg)

    C1("M <- (id - exponent * d)")
    M <- id - exponent * d ; C2()

    C1("b <- degree(graph,.)")
    b <- degree(graph, mode="out") ; C2()

    if(trace >= 2) {
	cat("  M : class ", class(M),";	   dim: ", dim(M),"\n")
	cat("  b : class ", class(b),"; length: ", length(b),"\n")
    }

    ## solve it

    ## MM: This is "the horror"	 -- ("it's the economy, stupid !!")
    ##	   particularly as  solve(M) is *never* sparse !!!
    ## He should use  solve(M, b) !!
    if(FALSE) {
	C1("solve(M) %*% b")
	ev <- solve(M, tol=tol) %*% b; C2()
    } else {
	C1("solve(M, b)")
	ev <- solve(M, b) ; C2()
    }

    if (rescale) {
	ev <- ev/sum(ev)
    } else {
	ev <- ev * sqrt(vg/sum((ev)^2))
    }

    ev[as.numeric(nodes) + 1]
}