[Rcpp-devel] Sparse matrix and RcppEigen

[Balamuta, James Joseph]

Greetings and Salutations Wagner,

1. I think you solved your own initial problem as it relates to solving matrices within eigen using the “solver” class.

e.g.

SparseMatrix<double> A;
SparseMatrix<double> B;

SolverClassName<SparseMatrix<double> > solver(A);

SparseMatrix<double> x = solver.solve(B);

2. Taking a square root of a matrix in Eigen is not supported officially. There does exist an unsupported square root feature for standard (e.g. not sparse matrices).

See: http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1MatrixSquareRoot.html

To use these features, add the include:  #include <unsupported/Eigen/MatrixFunctions>

3. Or write your own square root function using eigenvalues and eigenvectors!

Though, it seems that Eigen does not have built in support for obtaining eigenvalues & eigenvectors for a sparse matrix. You might need a secondary library to do so.

On the other hand, such support exists in Armadillo. However, ARPACK must be enabled / linked to. Unfortunately, this will probably not ever happen for RcppArmadillo since it is a pure template package.

With that being said, here’s some sample code for armadillo:

#include <RcppArmadillo.h>
using namespace Rcpp;

// Non-sparse decomp
// Assumes eigen decomp is real.
// [[Rcpp::export]]
arma::mat sqrt_mat(arma::mat x){

  arma::cx_vec eigval;
  arma::cx_mat eigvec;

  eig_gen(eigval, eigvec, x);

  return arma::conv_to<arma::mat>::from(
    eigvec*diagmat(1.0/sqrt(eigval))*trans(eigvec)
    );
}


// Sparse decomp
// Requires ARPACK to be enabled
// Assumes eigen decomp is real.
// [[Rcpp::export]]
arma::sp_mat sqrt_mat_sp(arma::sp_mat x){

  arma::cx_vec eigval;
  arma::cx_mat eigvec;

  eigs_gen(eigval, eigvec, x, x.n_cols);

  return
    eigvec*diagmat(1.0/sqrt(eigval))*trans(eigvec);
}

/*** R
x = matrix(c(5,2,2,5),nrow=2)

sqrt_mat(x)

library(Matrix)
sx = Matrix(X, sparse=T)

# Will error without ARPACK enabled.
#sqrt_mat_sp(x)
*/

Sincerely,

JJB


On Jun 8, 2015, at 2:37 AM, Wagner Bonat <wbonat at gmail.com<mailto:wbonat at gmail.com>> wrote:

Dear all,

I need to solve a linear system using only the lower triangular matrix from the Cholesky decomposition of a SPD matrix. My goal is to compute C^{-1/2} V^{-1/2} where V^{-1/2} and C are known matrices.
For general matrix class MatrixXd I can do using the code below. But, now I would like to do the same using the SparseMatrix class. I tried to adapt my code, but it did not work. Any help will be welcome.


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);'


srcsparse <- '
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::Lower> solver;
SparseMatrix<double>  LT = solver.compute(Omega).matrixL().solve(invSqrtV);
return wrap(LT);'

The problem is that there is no .solve() method associated with the matrixL().
It is just a intermediary step when computing

solver.compute(Omega).solve(invSqrtV)

This code computing C^{-1}V^{-1/2}, but I want only C^{-1/2}V^{-1/2}. It will be used on the Generalized Kronecker product in the context of multivariate spatial models.

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 at lists.r-forge.r-project.org<mailto:Rcpp-devel at lists.r-forge.r-project.org>
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

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