[Rcpp-devel] symmatu does not give symmetric matrix

[Jendoubi, Takoua]

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:bates at stat.wisc.edu]
Sent: 27 July 2016 14:45
To: Dirk Eddelbuettel; Jendoubi, Takoua
Cc: rcpp-devel at 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 <edd at debian.org<mailto:edd at 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 | edd at debian.org<mailto:edd at debian.org>
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel at lists.r-forge.r-project.org<mailto:Rcpp-devel at lists.r-forge.r-project.org>
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20160729/6c1eb6bb/attachment.html>