[Rcpp-devel] R & RcppArmadillo decomposition disagreement

[Dirk Eddelbuettel]

On 30 May 2013 at 20:46, Colin Rundel wrote:
| 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.

Not necessarily. You need to follow the actual codepath from R and Arma. 
 
| 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).

Well that can easily be the issue, and in that case you should probably
follow up with respective BLAS maintainers.  

The very nice things is that you plug BLAS in and out at essentially zero
cost (apart from the download and apt-get mechanics, say 15 secs each).  So
try

   reference blas (packages libblas3 and liblapack)
   atlas          (package libatlas3-base or a tuned variant)
   open-blas      (package libopenblas-base)  

See if one or more of these behave.

Dirk
 
| 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.
| 
| xuntyped binary data, test.R         [Click mouse-2 to save to a file]
| xuntyped binary data, util.cpp       [Click mouse-2 to save to a file]
| 
| 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/
| 
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com