<div dir="ltr"><div><div>Hello ! <br><br>Thank you for reply! <br><br></div>I tried to use your suggestions. The code compiled, but the result does not make sense. Have a look,<br><br></div>My really simple R function <br><div><p class="MsoNormal"><br><span style="font-size:11pt;font-family:"Calibri",sans-serif;color:rgb(31,73,125)"></span></p><p class="MsoNormal"><span style="font-size:11pt;font-family:"Calibri",sans-serif;color:rgb(31,73,125)">my.partialsolveR <- function(Omega, invSqrtV){solve(chol(Omega), invSqrtV)}</span></p><p class="MsoNormal"><br><span style="font-size:11pt;font-family:"Calibri",sans-serif;color:rgb(31,73,125)"></span></p><p class="MsoNormal">The C++ function using your advices</p><p class="MsoNormal"><br></p><p class="MsoNormal">partialsolveCS <- '<br>using namespace Rcpp;<br>using namespace Eigen;<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(Omega);<br>SparseMatrix<double> L = solver.matrixL();<br>Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver2(L);<br>SparseMatrix<double> L2 = solver2.solve(invSqrtV);<br>return wrap(L2);'<br><br>my.partialsolveCS <- cxxfunction(signature(Omegas = "dgCMatrix", invSqrtVs = "dgCMatrix"), <br></p><p class="MsoNormal"> body = partialsolveCS, plugin = "RcppEigen")<br></p><p class="MsoNormal"><br></p><p class="MsoNormal">I know that this type of conversion is not recommended but consider only for this example. <br></p><p class="MsoNormal"><br></p><p class="MsoNormal">Omega <- Matrix(c(1,0.8,0.8,1),2,2)<br>invSqrtV <- Diagonal(2,1)<br>Omega.M <- as(as.matrix(Omega),"dgCMatrix")<br>invSqrtV.M <- as(as.matrix(invSqrtV),"dgCMatrix")<br></p><p class="MsoNormal"><br></p><p class="MsoNormal">my.partialsolveR(Omega,invSqrtV)<br>2 x 2 Matrix of class "dgeMatrix"<br> [,1] [,2]<br>[1,] 1 -1.333333<br>[2,] 0 1.666667</p><p class="MsoNormal"><br></p><p class="MsoNormal">my.partialsolveCS(Omega.M, invSqrtV.M)<br>2 x 2 sparse Matrix of class "dgCMatrix"<br> <br>[1,] . .<br>[2,] 2.121996e-314 .<br></p><p class="MsoNormal"><br></p><p class="MsoNormal">As you can see there are something wrong in the code. I tried change Eigen::Lower for Eigen::Upper and vice-versa, but did not work.</p><p class="MsoNormal"><br></p><p class="MsoNormal">Thank you again !!<br> </p><br><br><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-06-10 5:48 GMT+02:00 Balamuta, James Joseph <span dir="ltr"><<a href="mailto:balamut2@illinois.edu" target="_blank">balamut2@illinois.edu</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div link="blue" vlink="purple" lang="EN-US">
<div>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">Regarding the code not compiling….<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">Try:<u></u><u></u></span></p><span class="">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">src <- '<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">using Eigen::Map;
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">using Eigen::SparseMatrix;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">using Eigen::LLT;<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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));<u></u><u></u></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));<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
</span><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)<u></u><u></u></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);<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">SparseMatrix<double> L = solver.matrixL();<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">// Second solver using the lower matrix<u></u><u></u></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);<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">SparseMatrix<double> L2 = solver2.solve(invSqrtV);<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">return wrap(L2);'<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">e.g.<u></u><u></u></span></p><span class="">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">Omega.m <- as.matrix(Omega)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">Vinvsqrt.m <- as.matrix(Vinvsqrt)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
</span><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.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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>
<u></u><u></u></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")<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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.”<u></u><u></u></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>
<u></u><u></u></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.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">Sincerely,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d">JJB<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri",sans-serif;color:#1f497d"><u></u> <u></u></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"> <a href="mailto:rcpp-devel-bounces@lists.r-forge.r-project.org" target="_blank">rcpp-devel-bounces@lists.r-forge.r-project.org</a> [mailto:<a href="mailto:rcpp-devel-bounces@lists.r-forge.r-project.org" target="_blank">rcpp-devel-bounces@lists.r-forge.r-project.org</a>]
<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> <a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank">rcpp-devel@lists.r-forge.r-project.org</a><br>
<b>Subject:</b> Re: [Rcpp-devel] Sparse matrix and RcppEigen<u></u><u></u></span></p><div><div class="h5">
<p class="MsoNormal"><u></u> <u></u></p>
<div>
<div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Dear all,<u></u><u></u></p>
</div>
<p class="MsoNormal">As requested I will try provide a better explanation. Consider the code below:<u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></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}.<u></u><u></u></p>
</div>
<p class="MsoNormal">Some auxiliary functions to provide an example.<u></u><u></u></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)<u></u><u></u></p>
</div>
<p class="MsoNormal"><u></u> <u></u></p>
<div>
<p class="MsoNormal">## R function<u></u><u></u></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")<u></u><u></u></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<u></u><u></u></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")<u></u><u></u></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<u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">Now, the results are completely different because my C++ code does not take into account the sparse structure.<u></u><u></u></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);'<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Any help is welcome. <u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">Thank you!<u></u><u></u></p>
</div>
</div>
</div>
</div>
</div></div><div>
<p class="MsoNormal"><br>
-- <u></u><u></u></p><span class="">
<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)<u></u><u></u></p>
</div>
</div>
</div>
</div>
</div>
</span></div>
</div>
</div>
<br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">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><br></blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>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)<br><br></div></div></div></div></div>
</div>