On 15/10/2013 20:37, Simon Urbanek 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?
That was my guess, and the workaround is in 3.0.2 so please try current
R (as the posting guide asked to be done before posting).
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
--
Brian D. Ripley, [email protected]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
_______________________________________________
R-SIG-Mac mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
