[Rcpp-devel] Sparse matrix and RcppEigen

[Balamuta, James Joseph]

Thanks for the reproducible example!

Here is the solution:

library(Matrix)

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.vector(diag(invSqrtV))
partialsolveCS <- '
using namespace Rcpp;
using namespace Eigen;
const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));
const VectorXd invSqrtV(as<VectorXd>(invSqrtVs));
MatrixXd mat = invSqrtV.asDiagonal();

Eigen::SimplicialLLT<SparseMatrix<double> > solver;

MatrixXd L = solver.compute(Omega).matrixL().solve(mat);

return wrap(L.sparseView());'

my.partialsolveCS <- cxxfunction(signature(Omegas = "dgCMatrix", invSqrtVs = "numeric"),
                                 body = partialsolveCS, plugin = "RcppEigen")


my.partial.solveR <- function(Omega,invSqrtV){ solve(t(chol(Omega)), invSqrtV)}

my_partialsolveCS(Omega.M, invSqrtV.M)
 my.partial.solveR(Omega, invSqrtV)


Few notes:

1.       The previous solution suffers from taking the cholesky decomp twice. Once on omega and then again on L matrix from Omega.  The solver.info() returned false since it could not do it. Anything further on was then zero’d out.

2.       The new approach uses the requirement that after the sparse solver has broken apart the matrix, it can act upon a matrix / vector that is dense.
(e.g. Ax = b => x = A^(-1)b, where A is sparse and b is dense)

3.       We will exploit the property that the second matrix is diagonal based by passing in a vector instead of an actual matrix and gaining the solve benefit of having the matrix registered as a diagonal matrix within eigen.

4.       The unfortunate side effect is the return of a matrix that is dense. We can easily switch between types using:

a.       SparseMatrix => MatrixXd: MatrixXd changed_to_dense = MatrixXd(sparse_matrix_to_change)

b.      MatrixXd => SparseMatrix: SparseMatrix<double> changed_to_sparse = dense_matrix_to_change.sparseView();

5.       We shouldn’t lose lots of speed using this approach due to how Eigen library is built (e.g. http://eigen.tuxfamily.org/dox/TopicInsideEigenExample.html )



Misc:
I’d highly encourage you in the future to use the sourceCpp form as it allows for comments and easy of reuse.

Rcpp::sourceCpp('path/to/file/eigen_sparse_solve.cpp')


#include <RcppEigen.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppEigen)]]

using Eigen::Map;
using Eigen::SparseMatrix;
using Eigen::LLT;

// [[Rcpp::export]]
Eigen::MatrixXd my_partialsolveCS( Eigen::SparseMatrix<double> Omega,  Eigen::VectorXd invSqrtV){

  // If we know invSqrtV is a diagonal, we'll use this to our advantage
  Eigen::MatrixXd mat = invSqrtV.asDiagonal();

  // Set up sparse solver.
  // Note, sparse solvers require one matrix to be sparse (omega) and the other matrix or vector to be dense.
  Eigen::SimplicialLLT<SparseMatrix<double> > solver;

  // Access the Lower Cholemsky result and use it to solve the system.
  Eigen::MatrixXd L = solver.compute(Omega).matrixL().solve(mat);

  // Convert to sparse matrix and return.
  return L.sparseView();
}

/*** R
# Runs the test code on compile
library(Matrix)

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.vector(diag(invSqrtV))


my.partial.solveR <- function(Omega,invSqrtV){ solve(t(chol(Omega)), invSqrtV)}

my.partial.solveR(Omega, invSqrtV)

my_partialsolveCS(Omega.M, invSqrtV.M)
*/


Sincerely,

JJB

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150611/bd060ebd/attachment-0001.html>