Hello Wagner,
I think there is some confusion in your question. From your first code
chunk, it seems that you first compute the Cholesky decomposition C=LL',
and then calculate L^{-1} * V^{-1/2}. However, this is *NOT* equal to C^{-1/2}
V^{-1/2}.
In the usual sense, C^{1/2} is defined to be a matrix M such that MM=C,
i.e., no transpose here. However, the Cholesky decomposition calculates
matrix L that satisfies LL'=C.
If C is symmetric, M is also symmetric but L is lower triangular.
Could you clarify your intention here?
Best,
Yixuan
2015-06-08 3:37 GMT-04:00 Wagner Bonat <wbo...@gmail.com>:
> 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
>
--
Yixuan Qiu <yixuan....@cos.name>
Department of Statistics,
Purdue University
_______________________________________________
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
