<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoPlainText, li.MsoPlainText, div.MsoPlainText
{mso-style-priority:99;
mso-style-link:"Plain Text Char";
margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
span.EmailStyle17
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.PlainTextChar
{mso-style-name:"Plain Text Char";
mso-style-priority:99;
mso-style-link:"Plain Text";
font-family:"Calibri","sans-serif";}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:72.0pt 72.0pt 72.0pt 72.0pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-GB" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoPlainText">Dear all,<o:p></o:p></p>
<p class="MsoPlainText"><o:p> </o:p></p>
<p class="MsoPlainText">Thank you for your help.<o:p></o:p></p>
<p class="MsoPlainText">I am actually dealing with high dimensional data and using cholesky decomposition to compute the determinant of the matrix.<o:p></o:p></p>
<p class="MsoPlainText">The example I gave was only intended to be illustrative and it brought some confusion.<o:p></o:p></p>
<p class="MsoPlainText"><o:p> </o:p></p>
<p class="MsoPlainText">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.<o:p></o:p></p>
<p class="MsoPlainText"><o:p> </o:p></p>
<p class="MsoPlainText">Thank you again,<o:p></o:p></p>
<p class="MsoPlainText">Takoua.<o:p></o:p></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><b><span lang="EN-US" style="font-size:10.0pt;font-family:"Tahoma","sans-serif"">From:</span></b><span lang="EN-US" style="font-size:10.0pt;font-family:"Tahoma","sans-serif""> Douglas Bates [mailto:bates@stat.wisc.edu]
<br>
<b>Sent:</b> 27 July 2016 14:45<br>
<b>To:</b> Dirk Eddelbuettel; Jendoubi, Takoua<br>
<b>Cc:</b> rcpp-devel@lists.r-forge.r-project.org; Martin Maechler<br>
<b>Subject:</b> Re: [Rcpp-devel] symmatu does not give symmetric matrix<o:p></o:p></span></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal">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.<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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 <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">prod(x)<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">The dput function will print the values out to full precision so try<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">dput(prod(x))<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">dput(x %*% t(x))<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">to see the full precision results.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">dput(tcrossprod(x))<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">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.<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal">On Wed, Jul 27, 2016 at 7:40 AM Dirk Eddelbuettel <<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a>> wrote:<o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0cm 0cm 0cm 6.0pt;margin-left:4.8pt;margin-right:0cm">
<p class="MsoNormal"><br>
Sorry for the delay in responding. Mail delivery was stuck at R-Forge for a<br>
few days it seems.<br>
<br>
On 21 July 2016 at 10:46, Jendoubi, Takoua wrote:<br>
| Dear all,<br>
|<br>
| I am using RcppArmadillo to deal with some matrix computations. Specifically I<br>
| need to cholesky factorization of some symmetric matrices.<br>
|<br>
| I am generating random vectors using Rcpp and using them to construct symmetric<br>
| matrices.<br>
|<br>
| I have an error stating that my matrix is not symmetric although it definitely<br>
| should be. Here is the example I am working with (in R):<br>
|<br>
| >x<br>
| [1] -1.6683320 -0.8597148<br>
| >x%*%t(x)<br>
| [,1] [,2]<br>
| [1,] 2.783332 1.4342896<br>
| [2,] 1.434290 0.7391095<br>
|<br>
| Apparently, it is a rounding-off error. Is there any way to ensure that x%*%t<br>
| (x) gives an exactly symmetric matrix to use for cholesky factorization?<br>
<br>
I don't know of an automatic way. R surely has no native type for symmetric<br>
matrices (at least not at the SEXP level as it really only has vectors with<br>
dim attributes). I would probably start with some helper functions at the R<br>
or C++ level. Some may have better tricks...<br>
<br>
Dirk<br>
<br>
<br>
--<br>
<a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a> | @eddelbuettel |
<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a><br>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><o:p></o:p></p>
</blockquote>
</div>
</div>
</div>
</body>
</html>