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