[Rcpp-devel] Sparse matrix and RcppEigen

[Balamuta, James Joseph]

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-bounces at lists.r-forge.r-project.org [mailto:rcpp-devel-bounces at 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 at 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)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150610/3662e367/attachment-0001.html>