<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 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
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.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0in;
margin-right:0in;
margin-bottom:0in;
margin-left:.5in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman",serif;}
span.hoenzb
{mso-style-name:hoenzb;}
span.EmailStyle19
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:201797004;
mso-list-type:hybrid;
mso-list-template-ids:-1472038214 67698703 67698713 67698715 67698703 67698713 67698715 67698703 67698713 67698715;}
@list l0:level1
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></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-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Regarding the code not compiling….<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Try:<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">src <- '<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">using Eigen::Map;
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">using Eigen::SparseMatrix;<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">using Eigen::LLT;<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">// First solver associated with Omega matrix (taken w.r.t to Upper triangular view)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Eigen::SimplicialLLT<SparseMatrix<double> , Eigen::Upper > solver(Omega);<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">SparseMatrix<double> L = solver.matrixL();<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">// Second solver using the lower matrix<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Eigen::SimplicialLLT<SparseMatrix<double> , Eigen::Lower > solver2(L);<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">SparseMatrix<double> L2 = solver2.solve(invSqrtV);<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">return wrap(L2);'<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"><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">The other issue is how the matrices are passed in.<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Before passing into C++, you are converting the objects into actual R matrices without any sparse properties:<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">e.g.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Omega.m <- as.matrix(Omega)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Vinvsqrt.m <- as.matrix(Vinvsqrt)<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">The matrices must be of the appropriate sparse type in R to be translated into sparse C++ objects.<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">This leads to the slots being formatted as:<br>
<br>
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">my.partial.solveC <- cxxfunction(signature(Omegas = "dgCMatrix", invSqrtVs = " ddiMatrix "), body = src, plugin = "RcppEigen")<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"><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">However, there is an issue with how the second matrix is typed.<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Specifically, it is typed as: ddiMatrix.
<br>
<br>
This leads to an error when calling the compile statement of: “Error: No such slot.”<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D"><br>
So, I have a funny feeling the conversions in place do not account for the sparse diagonal matrix class (ddiMatrix).<br>
<br>
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">It is at this point, that I can’t advise you anymore outside of saying you may want to consider building the invSqrtV and Omega matrices within Eigen to avoid
this conversion issue.<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">Sincerely,<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D">JJB<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"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri",sans-serif">From:</span></b><span style="font-size:11.0pt;font-family:"Calibri",sans-serif"> rcpp-devel-bounces@lists.r-forge.r-project.org [mailto:rcpp-devel-bounces@lists.r-forge.r-project.org]
<b>On Behalf Of </b>Wagner Bonat<br>
<b>Sent:</b> Tuesday, June 09, 2015 6:50 AM<br>
<b>To:</b> Yixuan Qiu<br>
<b>Cc:</b> rcpp-devel@lists.r-forge.r-project.org<br>
<b>Subject:</b> Re: [Rcpp-devel] Sparse matrix and RcppEigen<o:p></o:p></span></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Dear all,<o:p></o:p></p>
</div>
<p class="MsoNormal">As requested I will try provide a better explanation. Consider the code below:<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Let C = V^{1/2}%*%Omega%*%V^{1/2} and C^{1/2} be the Cholesky decomposition, such that C^{1/2}%*%C^{1/2}^T = C.<br>
My goal is to compute Omega^{-1/2}%*%V{-1/2}. <br>
My idea is to use that I know V^{-1/2} and rewrite C^{-1} = V^{-1/2}%*%Omega^{-1}%*%V^{-1/2} = V^{-1/2}%*%Omega^{-1/2}%*%(Omega^{-1/2}%*%V^{-1/2})^T<br>
The main term now is Omega^{-1/2}%*%V^{-1/2} for computing it my idea is first compute the Cholesky denoted by Omega^{1/2} and then using the triangular structure to solve a system to get Omega^{-1/2}%*%V^{-1/2}.<o:p></o:p></p>
</div>
<p class="MsoNormal">Some auxiliary functions to provide an example.<o:p></o:p></p>
<div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt"><br>
# Inverse of square root of V<br>
build.invSqrtV <- function(mu, power,n.sample){Diagonal(n.sample, 1/sqrt(mu^power))}<br>
<br>
## Omega matrix - It is an example of dense structure.<br>
build.Omega <- function(par,U){par[1]*exp(-U/par[2])}<br>
<br>
## Small sample size<br>
x1 <- seq(0,1,l=30)<br>
x2 <- seq(0,1,l=30)<br>
grid <- expand.grid(x1,x2)<br>
U <- forceSymmetric(as.matrix(dist(grid)))<br>
cov <- seq(-1,1,l=900)<br>
X <- model.matrix(~ cov)<br>
mu <- exp(X%*%c(1,0.5))<br>
<br>
Vinvsqrt <- build.invSqrtV(mu, power = 2, n.sample =900)<br>
Omega <- build.Omega(par=c(1,0.25),U)<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal">## R function<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">my.partial.solveR <- function(Omega,invSqrtV){ solve(t(chol(Omega)), invSqrtV)}<br>
result1 <- my.partial.solveR(Omega, Vinvsqrt)<br>
<br>
## Now using RcppEigen<br>
src <- '<br>
using namespace Rcpp;<br>
using namespace Eigen;<br>
const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas));<br>
const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs));<br>
MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV); <br>
return wrap(LT);'<br>
<br>
## Here I need to convert from dsyMatrix to matrix<br>
Omega.m <- as.matrix(Omega)<br>
Vinvsqrt.m <- as.matrix(Vinvsqrt)<br>
my.partial.solveC <- cxxfunction(signature(Omegas = "mat", invSqrtVs = "mat"), body = src, plugin = "RcppEigen")<br>
result2 <- my.partial.solveC(Omega.m,Vinvsqrt.m)<br>
<br>
## Comparing the results<br>
tt <- abs(result1) - abs(result2)<br>
sum(tt)<br>
<br>
## Comparing the computational time<br>
benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m, Vinvsqrt.m), order = "relative")<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt"> test replications elapsed relative<br>
2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 1.000<br>
1 my.partial.solveR(Omega, Vinvsqrt) 100 25.414 2.855<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">You can see that the C++ code is much faster than my R code.<br>
<br>
But, now consider that my Omega matrix is sparse<br>
z <- rep(1,10)<br>
Z <- z%*%t(z)<br>
listZ <- list()<br>
for(i in 1:90){listZ[[i]] <- Z}<br>
Omega <- Diagonal(900,1) + 0.8*bdiag(listZ)<br>
Omega.m <- as.matrix(Omega)<br>
<br>
benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m, Vinvsqrt.m), order = "relative")<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><br>
test replications elapsed relative<br>
1 my.partial.solveR(Omega, Vinvsqrt) 100 0.101 1.000<br>
2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 88.129<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Now, the results are completely different because my C++ code does not take into account the sparse structure.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">My idea is reproduze the same idea using SparseMatrix class. I try the following function, but it does not work, since there is no solve method associated with the solver. So, I do not know how to use the triangular structure to solve the
linear system<br>
<br>
src <- '<br>
using Eigen::Map; <br>
using Eigen::SparseMatrix;<br>
using Eigen::LLT;<br>
<br>
const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));<br>
const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));<br>
Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Upper> solver;<br>
solver.compute(Omega);<br>
SparseMatrix<double> L = solver.matrixL().solve(invSqrtV); #### Do not work<br>
return wrap(L);'<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Any help is welcome. <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">Thank you!<o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<p class="MsoNormal"><br>
-- <o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Wagner Hugo Bonat<br>
----------------------------------------------------------------------------------------------<br>
Department of Mathematics and Computer Science (IMADA)<br>
University of Southern Denmark (SDU) and<br>
Laboratório de Estatística e Geoinformação (LEG)<br>
Universidade Federal do Paraná (UFPR)<o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>