[R-SIG-Mac] Erroneous zeros in results using Matrix package

[Braun, Michael]

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.

On Oct 15, 2013, at 1:26 PM, Roger Koenker <rkoenker at illinois.edu> 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    rkoenker at uiuc.edu            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
>> R-SIG-Mac at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mac
>