I was unaware of the 32-bit overflow issue, and I must have missed the announcement of R-3.0.2. I cannot update from 3.0.1 to 3.0.2 right away, for some unrelated reasons, but I will do that ASAP. I did peruse the NEWS file for R-3.0.2, R-3.0.2-patched, and R-devel, but I did not see any items that immediately jumped out as being connected.
My installation of R links to the Intel MKL BLAS, version 11.0, Update 3. I compiled R using clang, with flags -msse4.2 -mtune=corei7 -O2 (among some others). I would doubt that the BLAS choice is the problem, since a simple coercion of a large identity matrix to its sparse counterpart should not have to call a BLAS routine. I can send my Makeconf file if you think that would be useful. On Oct 15, 2013, at 2:37 PM, Simon Urbanek <[email protected]> wrote: > 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 >> > -------------------------- Michael Braun Associate Professor of Marketing Cox School of Business Southern Methodist University Dallas, TX 75275 [email protected] _______________________________________________ R-SIG-Mac mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-mac
