[Rcpp-devel] R & RcppArmadillo decomposition disagreement

[Colin Rundel]

Hi everyone,

I have recently run into a strange situation with RcppArmadillo and was hoping someone might have some insight. I am working with a moderately sized (820x820) covariance matrix in C++ which causes arms::chol to fail. In trying to diagnose the issue I have written out the offending matrix using the arma_binary format and then read it back into R using a small helper function. Once in R, I am able to use the R chol function to decompose the matrix without issue (no pivoting).

This makes no sense to me, as it is my understanding that both R and armadillo should be using the system lapack library and hence both should call the same dpotrf. To make matters worse, I've also tried the eigen decomposition of the same matrix, which results in anything from a crash / system hang / garbage when using armadillo (eig_sym) but the same matrix read into R will happily return a valid result from the eigen function. Other, usually smaller, matrices work fine on the same system.

One further wrinkle, is that this also seems to be somewhat system dependent as it occurs on my Ubuntu box (13.04, R 3.0.1, RcppArmadillo 0.3.820) but not my Macbook (10.8.3, R 3.0.1, RcppArmadillo 0.3.820). This leads me to believe that the issue may be with lapack/blas but I see the same result when using base lapack, atlas (3.8.4-9ubuntu1), or openblas (0.2.6-1~exp1).

I have attached a minimal working example that replicates the issue for my system, and I have also included a link to the matrix file as it is somewhat large.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.R
Type: application/octet-stream
Size: 218 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130530/ab518fc4/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: util.cpp
Type: application/octet-stream
Size: 1379 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130530/ab518fc4/attachment-0001.obj>
-------------- next part --------------

https://www.dropbox.com/s/6i1c8xp2dsthq4y/tmp.dat

-Colin

-----

Colin Rundel
Postdoctoral Associate
Duke University, Department of Statistical Science
colin.rundel at stat.duke.edu
www.stat.duke.edu/~cr173/