On Oct 15, 2013, at 3:10 PM, Braun, Michael wrote: > I suppose that is true. But the problem does not occur in smaller examples. > It only occurs when working with very large matrices. That's the whole > point. And to be fair, I did warn the members of the list that this is a > memory-intensive example. >
Isn't this related to the mysterious issue with 32-bit overflow in matrices on some OS X systems? Unfortunately I can't seem to find that thread - I recall Brian committing a work-around, so can you test more recent R? In addition, which BLAS implementation are you using? Thanks, Simon > On Oct 15, 2013, at 1:26 PM, Roger Koenker <[email protected]> wrote: > >> You would be much more likely to get a response to this if you had a >> _smaller_ reproducible example. >> >> url: www.econ.uiuc.edu/~roger Roger Koenker >> email [email protected] Department of Economics >> vox: 217-333-4558 University of Illinois >> fax: 217-244-6678 Urbana, IL 61801 >> >> On Oct 15, 2013, at 11:36 AM, Braun, Michael wrote: >> >>> I've come across some interesting behavior using the Matrix package, and I >>> cannot seem to explain it. I have sent this issue to the Matrix >>> maintainers, and Martin Maechler was unable to replicate it on his Linux >>> system. Since it might be something Mac-specific, he suggested that I post >>> here. The complete code is pasted at the bottom of this email, along with >>> my sessionInfo() output. >>> >>> The function f takes 3 parameters: p, k and N. First, it creates A, a >>> lower triangle of a sparse block diagonal matrix with p blocks, and each >>> block being k x k. I then let S = AA', and have L be the lower Cholesky >>> decomposition of L. (I use Cholesky instead of chol because, for my "real" >>> application, because I need the permutation matrix also). I then multiply >>> L by a matrix z, where each column in z is a standard MVN sample. >>> >>> None of the elements should be zero. However, if z becomes very large, I >>> get lots of zero entries. >>> >>> Additionally, you will see that when I try to coerce a large identity >>> matrix to a dtCMatrix, some of the ones are lost. >>> >>> A word of warning: this script does consume a lot of RAM, and it will take >>> several minutes to run. But it would be helpful if someone with a Mac Pro >>> with 32GB of RAM could try to replicate the problem. >>> >>> Thanks in advance for your help. >>> >>> Michael >>> >>> >>> library(Matrix) >>> >>> f <- function(p,k,N) { >>> A <- tril(kronecker(Diagonal(p),Matrix(1:(k^2),k,k))) >>> S <- tcrossprod(A) >>> CH <- Cholesky(S) >>> L <- expand(CH)$L >>> z <- rnorm(N*k*p) >>> dim(z) <- c(p*k,N) >>> y <- L %*% z >>> return(sum(y==0)) >>> } >>> >>> system.time(small <- f(p= 100, k=3, N= 200) ) >>> system.time(medium <- f(p= 200, k=3, N= 5000) ) >>> >>> print(small ) ## 0 >>> print(medium) ## 0 >>> >>> print( system.time(large.p <- f(p=50000, k=3, N= 5000) )) >>> print( system.time(large.N <- f(p=10000, k=3, N=20000) )) >>> >>> print(large.p) ## I get 536870912, but it should be 0 >>> print(large.N) ## I get 536870912, but it should be 0 >>> >>> system.time( W1 <- as(diag(10000),"dtCMatrix") ) >>> print(sum(diag(W1))) ## should be 10000, and that's what I get >>> >>> system.time( W3 <- as(diag(30000),"dtCMatrix") ) >>> print(sum(diag(W3)))## should be 30000, I get 12104 >>> >>> sessionInfo() >>> R version 3.0.1 (2013-05-16) >>> Platform: x86_64-apple-darwin12.3.0/x86_64 (64-bit) >>> >>> locale: >>> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 >>> >>> attached base packages: >>> [1] stats graphics grDevices utils datasets methods base >>> >>> other attached packages: >>> [1] Matrix_1.0-14 lattice_0.20-15 >>> >>> loaded via a namespace (and not attached): >>> [1] grid_3.0.1 tools_3.0.1 >>> >>> -------------------------- >>> Michael Braun >>> Associate Professor of Marketing >>> Cox School of Business >>> Southern Methodist University >>> Dallas, TX 75275 >>> braunm _at_ smu.edu >>> >>> _______________________________________________ >>> R-SIG-Mac mailing list >>> [email protected] >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mac >> > > _______________________________________________ > R-SIG-Mac mailing list > [email protected] > https://stat.ethz.ch/mailman/listinfo/r-sig-mac > _______________________________________________ R-SIG-Mac mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-mac
