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
_______________________________________________
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
