Hello !
Thank you for reply!
I tried to use your suggestions. The code compiled, but the result does not
make sense. Have a look,
My really simple R function
my.partialsolveR <- function(Omega, invSqrtV){solve(chol(Omega), invSqrtV)}
The C++ function using your advices
partialsolveCS <- '
using namespace Rcpp;
using namespace Eigen;
const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));
const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));
Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Upper> solver(Omega);
SparseMatrix<double> L = solver.matrixL();
Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver2(L);
SparseMatrix<double> L2 = solver2.solve(invSqrtV);
return wrap(L2);'
my.partialsolveCS <- cxxfunction(signature(Omegas = "dgCMatrix", invSqrtVs
= "dgCMatrix"),
body = partialsolveCS,
plugin = "RcppEigen")
I know that this type of conversion is not recommended but consider only
for this example.
Omega <- Matrix(c(1,0.8,0.8,1),2,2)
invSqrtV <- Diagonal(2,1)
Omega.M <- as(as.matrix(Omega),"dgCMatrix")
invSqrtV.M <- as(as.matrix(invSqrtV),"dgCMatrix")
my.partialsolveR(Omega,invSqrtV)
2 x 2 Matrix of class "dgeMatrix"
[,1] [,2]
[1,] 1 -1.333333
[2,] 0 1.666667
my.partialsolveCS(Omega.M, invSqrtV.M)
2 x 2 sparse Matrix of class "dgCMatrix"
[1,] . .
[2,] 2.121996e-314 .
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.
Thank you again !!
2015-06-10 5:48 GMT+02:00 Balamuta, James Joseph <balam...@illinois.edu>:
> 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
>
--
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
