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