Regarding the code not compiling…. Try:
src <- ' using Eigen::Map; using Eigen::SparseMatrix; using Eigen::LLT; const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas)); const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs)); // First solver associated with Omega matrix (taken w.r.t to Upper triangular view) Eigen::SimplicialLLT<SparseMatrix<double> , Eigen::Upper > solver(Omega); SparseMatrix<double> L = solver.matrixL(); // Second solver using the lower matrix Eigen::SimplicialLLT<SparseMatrix<double> , Eigen::Lower > solver2(L); SparseMatrix<double> L2 = solver2.solve(invSqrtV); return wrap(L2);' The other issue is how the matrices are passed in. Before passing into C++, you are converting the objects into actual R matrices without any sparse properties: e.g. Omega.m <- as.matrix(Omega) Vinvsqrt.m <- as.matrix(Vinvsqrt) The matrices must be of the appropriate sparse type in R to be translated into sparse C++ objects. This leads to the slots being formatted as: my.partial.solveC <- cxxfunction(signature(Omegas = "dgCMatrix", invSqrtVs = " ddiMatrix "), body = src, plugin = "RcppEigen") However, there is an issue with how the second matrix is typed. Specifically, it is typed as: ddiMatrix. This leads to an error when calling the compile statement of: “Error: No such slot.” So, I have a funny feeling the conversions in place do not account for the sparse diagonal matrix class (ddiMatrix). 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. Sincerely, JJB From: rcpp-devel-boun...@lists.r-forge.r-project.org [mailto:rcpp-devel-boun...@lists.r-forge.r-project.org] On Behalf Of Wagner Bonat Sent: Tuesday, June 09, 2015 6:50 AM To: Yixuan Qiu Cc: rcpp-devel@lists.r-forge.r-project.org Subject: Re: [Rcpp-devel] Sparse matrix and RcppEigen Dear all, As requested I will try provide a better explanation. Consider the code below: 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. My goal is to compute Omega^{-1/2}%*%V{-1/2}. 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 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}. Some auxiliary functions to provide an example. # Inverse of square root of V build.invSqrtV <- function(mu, power,n.sample){Diagonal(n.sample, 1/sqrt(mu^power))} ## Omega matrix - It is an example of dense structure. build.Omega <- function(par,U){par[1]*exp(-U/par[2])} ## Small sample size x1 <- seq(0,1,l=30) x2 <- seq(0,1,l=30) grid <- expand.grid(x1,x2) U <- forceSymmetric(as.matrix(dist(grid))) cov <- seq(-1,1,l=900) X <- model.matrix(~ cov) mu <- exp(X%*%c(1,0.5)) Vinvsqrt <- build.invSqrtV(mu, power = 2, n.sample =900) Omega <- build.Omega(par=c(1,0.25),U) ## R function my.partial.solveR <- function(Omega,invSqrtV){ solve(t(chol(Omega)), invSqrtV)} result1 <- my.partial.solveR(Omega, Vinvsqrt) ## Now using RcppEigen src <- ' using namespace Rcpp; using namespace Eigen; const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas)); const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs)); MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV); return wrap(LT);' ## Here I need to convert from dsyMatrix to matrix Omega.m <- as.matrix(Omega) Vinvsqrt.m <- as.matrix(Vinvsqrt) my.partial.solveC <- cxxfunction(signature(Omegas = "mat", invSqrtVs = "mat"), body = src, plugin = "RcppEigen") result2 <- my.partial.solveC(Omega.m,Vinvsqrt.m) ## Comparing the results tt <- abs(result1) - abs(result2) sum(tt) ## Comparing the computational time benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m, Vinvsqrt.m), order = "relative") test replications elapsed relative 2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 1.000 1 my.partial.solveR(Omega, Vinvsqrt) 100 25.414 2.855 You can see that the C++ code is much faster than my R code. But, now consider that my Omega matrix is sparse z <- rep(1,10) Z <- z%*%t(z) listZ <- list() for(i in 1:90){listZ[[i]] <- Z} Omega <- Diagonal(900,1) + 0.8*bdiag(listZ) Omega.m <- as.matrix(Omega) benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m, Vinvsqrt.m), order = "relative") test replications elapsed relative 1 my.partial.solveR(Omega, Vinvsqrt) 100 0.101 1.000 2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 88.129 Now, the results are completely different because my C++ code does not take into account the sparse structure. 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 src <- ' using Eigen::Map; using Eigen::SparseMatrix; using Eigen::LLT; const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas)); const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs)); Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Upper> solver; solver.compute(Omega); SparseMatrix<double> L = solver.matrixL().solve(invSqrtV); #### Do not work return wrap(L);' Any help is welcome. Thank you! -- Wagner Hugo Bonat ---------------------------------------------------------------------------------------------- Department of Mathematics and Computer Science (IMADA) University of Southern Denmark (SDU) and Laboratório de Estatística e Geoinformação (LEG) Universidade Federal do Paraná (UFPR)
_______________________________________________ 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
