Dear all,
Thank you for your help. I am actually dealing with high dimensional data and using cholesky decomposition to compute the determinant of the matrix. The example I gave was only intended to be illustrative and it brought some confusion. It turned out that the error (“matrix is not symmetric positive definite” ) is raised because the matrix became ill-conditioned during sampling and was containing very small values. The problem is solved by regularizing the matrix. Thank you again, Takoua. From: Douglas Bates [mailto:ba...@stat.wisc.edu] Sent: 27 July 2016 14:45 To: Dirk Eddelbuettel; Jendoubi, Takoua Cc: rcpp-devel@lists.r-forge.r-project.org; Martin Maechler Subject: Re: [Rcpp-devel] symmatu does not give symmetric matrix First, this is not a good way of constructing a matrix to test the Cholesky decomposition because the result is rank-deficient. It is a 2 by 2 matrix of rank 1. In fact it is the very definition of a rank-1 matrix - a column vector multiplied by its transpose. I am surprised at the result because I can't think of a way in computer arithmetic that the off-diagonals would be different. They should both be prod(x) The dput function will print the values out to full precision so try dput(prod(x)) dput(x %*% t(x)) to see the full precision results. Thirdly, R has two functions, crossprod and tcrossprod, to compute X'X and XX', respectively. They compute one triangle of the result and copy it into the other, so the result is guaranteed to be symmetric. Try dput(tcrossprod(x)) Finally, what are you going to do with the Cholesky decomposition? If you are doing this simply to get the decomposition, compiling C++ code, even with all the help provided by Rcpp, is doing things the hard way. Just use the chol function in R or the cholesky decomposition in the Matrix package if you want more control. On Wed, Jul 27, 2016 at 7:40 AM Dirk Eddelbuettel <e...@debian.org<mailto:e...@debian.org>> wrote: Sorry for the delay in responding. Mail delivery was stuck at R-Forge for a few days it seems. On 21 July 2016 at 10:46, Jendoubi, Takoua wrote: | Dear all, | | I am using RcppArmadillo to deal with some matrix computations. Specifically I | need to cholesky factorization of some symmetric matrices. | | I am generating random vectors using Rcpp and using them to construct symmetric | matrices. | | I have an error stating that my matrix is not symmetric although it definitely | should be. Here is the example I am working with (in R): | | >x | [1] -1.6683320 -0.8597148 | >x%*%t(x) | [,1] [,2] | [1,] 2.783332 1.4342896 | [2,] 1.434290 0.7391095 | | Apparently, it is a rounding-off error. Is there any way to ensure that x%*%t | (x) gives an exactly symmetric matrix to use for cholesky factorization? I don't know of an automatic way. R surely has no native type for symmetric matrices (at least not at the SEXP level as it really only has vectors with dim attributes). I would probably start with some helper functions at the R or C++ level. Some may have better tricks... Dirk -- http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org<mailto:e...@debian.org> _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org<mailto:Rcpp-devel@lists.r-forge.r-project.org> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
