<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">
Greetings and Salutations Wagner,
<div class=""><br class="">
</div>
<div class="">1. I think you solved your own initial problem as it relates to solving matrices within eigen using the “solver” class. </div>
<div class=""><br class="">
</div>
<div class="">e.g.</div>
<div class=""><br class="">
</div>
<div class="">SparseMatrix<double> A;</div>
<div class="">SparseMatrix<double> B;</div>
<div class=""><br class="">
</div>
<span class="">SolverClassName<SparseMatrix<double> > solver(A);<br class="">
</span><span class=""><br class="">
</span>SparseMatrix<double> <span class="">x = solver.solve(B);</span>
<div class=""><span class=""><br class="">
</span>
<div class="">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). </div>
<div class=""><br class="">
</div>
<div class="">See: <a href="http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1MatrixSquareRoot.html" class="">http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1MatrixSquareRoot.html</a></div>
<div class=""><br class="">
</div>
<span class="">To use these features, add the include: #include <unsupported/Eigen/MatrixFunctions></span></div>
<div class=""><span class=""><br class="">
</span></div>
<div class="">3. Or write your own square root function using eigenvalues and eigenvectors!</div>
<div class=""><br class="">
</div>
<div class="">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. </div>
<div class=""><br class="">
</div>
<div class="">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. </div>
<div class=""><br class="">
</div>
<div class="">With that being said, here’s some sample code for armadillo:</div>
<div class=""><br class="">
</div>
<div class="">
<div class="">#include <RcppArmadillo.h></div>
<div class="">using namespace Rcpp;</div>
<div class=""><br class="">
</div>
<div class="">// Non-sparse decomp</div>
<div class="">// Assumes eigen decomp is real.</div>
<div class="">// [[Rcpp::export]]</div>
<div class="">arma::mat sqrt_mat(arma::mat x){</div>
<div class=""> </div>
<div class=""> arma::cx_vec eigval;</div>
<div class=""> arma::cx_mat eigvec;</div>
<div class=""> </div>
<div class=""> eig_gen(eigval, eigvec, x);</div>
<div class=""> </div>
<div class=""> return arma::conv_to<arma::mat>::from(</div>
<div class=""> eigvec*diagmat(1.0/sqrt(eigval))*trans(eigvec)</div>
<div class=""> );</div>
<div class="">}</div>
</div>
<div class=""><br class="">
</div>
<div class=""><br class="">
</div>
<div class="">
<div class="">
<div class="">// Sparse decomp</div>
<div class="">// Requires ARPACK to be enabled</div>
<div class="">// Assumes eigen decomp is real.</div>
<div class="">// [[Rcpp::export]]</div>
<div class="">arma::sp_mat sqrt_mat_sp(arma::sp_mat x){</div>
<div class=""> </div>
<div class=""> arma::cx_vec eigval;</div>
<div class=""> arma::cx_mat eigvec;</div>
<div class=""> </div>
<div class=""> eigs_gen(eigval, eigvec, x, x.n_cols);</div>
<div class=""> </div>
<div class=""> return </div>
<div class=""> eigvec*diagmat(1.0/sqrt(eigval))*trans(eigvec);</div>
<div class="">}</div>
</div>
</div>
<div class=""><br class="">
</div>
<div class="">/*** R</div>
<div class="">
<div class="">x = matrix(c(5,2,2,5),nrow=2) </div>
<div class=""><br class="">
</div>
<div class="">sqrt_mat(x)</div>
<div class=""><br class="">
</div>
<div class="">library(Matrix)</div>
<div class="">sx = Matrix(X, sparse=T)</div>
</div>
<div class=""><br class="">
</div>
<div class=""># Will error without ARPACK enabled.</div>
<div class="">
<div class="">#sqrt_mat_sp(x)</div>
</div>
<div class="">*/</div>
<div class=""><br class="">
</div>
<div class="">
<div class="">Sincerely,</div>
<div class=""><br class="">
</div>
<div class="">JJB</div>
<div class="">
<div class=""><br class="">
<div>
<blockquote type="cite" class="">
<div class="">On Jun 8, 2015, at 2:37 AM, Wagner Bonat <<a href="mailto:wbonat@gmail.com" class="">wbonat@gmail.com</a>> wrote:</div>
<br class="Apple-interchange-newline">
<div class="">
<div dir="ltr" class="">
<div class="">Dear all,<br class="">
<br class="">
</div>
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.
<br class="">
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.
<div class="">
<div class=""><br class="">
<br class="">
src <- '<br class="">
using namespace Rcpp;<br class="">
using namespace Eigen;<br class="">
const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas));<br class="">
const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs));<br class="">
MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV);<br class="">
return wrap(LT);'<br class="">
<br class="">
<br class="">
srcsparse <- '<br class="">
using Eigen::Map; <br class="">
using Eigen::SparseMatrix;<br class="">
using Eigen::LLT;<br class="">
<br class="">
const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));<br class="">
const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));<br class="">
Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver;<br class="">
SparseMatrix<double> LT = solver.compute(Omega).matrixL()<b class="">.solve(invSqrtV);</b><br class="">
return wrap(LT);'<br class="">
<br class="">
</div>
<div class="">The problem is that there is no .solve() method associated with the matrixL().
<br class="">
</div>
<div class="">It is just a intermediary step when computing<br class="">
<br class="">
solver.compute(Omega).solve(invSqrtV)<br class="">
<br class="">
</div>
<div class="">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.<br class="">
<br class="">
</div>
<div class="">Thank you.<br class="">
</div>
<div class="">-- <br class="">
<div class="gmail_signature">
<div dir="ltr" class="">
<div class="">
<div dir="ltr" class="">
<div class="">Wagner Hugo Bonat<br class="">
----------------------------------------------------------------------------------------------<br class="">
Department of Mathematics and Computer Science (IMADA)<br class="">
University of Southern Denmark (SDU) and<br class="">
Laboratório de Estatística e Geoinformação (LEG)<br class="">
Universidade Federal do Paraná (UFPR)<br class="">
<br class="">
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
_______________________________________________<br class="">
Rcpp-devel mailing list<br class="">
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" class="">Rcpp-devel@lists.r-forge.r-project.org</a><br class="">
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</div>
</blockquote>
</div>
<br class="">
</div>
</div>
</div>
</body>
</html>