<div dir="ltr">Hello Wagner,<div>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} * <span style="font-size:14px">V^{-1/2}. However, this is <b>NOT</b> equal to </span><span style="font-size:14px">C^{-1/2} V^{-1/2}.</span></div><div><span style="font-size:14px">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.</span></div><div><span style="font-size:14px">If C is symmetric, M is also symmetric but L is lower triangular.</span></div><div><span style="font-size:14px"><br></span></div><div><span style="font-size:14px">Could you clarify your intention here?</span></div><div><span style="font-size:14px"><br></span></div><div><span style="font-size:14px"><br></span></div><div><span style="font-size:14px">Best,</span></div><div><span style="font-size:14px">Yixuan</span></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-06-08 3:37 GMT-04:00 Wagner Bonat <span dir="ltr"><<a href="mailto:wbonat@gmail.com" target="_blank">wbonat@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Dear all,<br><br></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>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><div><br><br>src <- '<br>using namespace Rcpp;<br>using namespace Eigen;<br>const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas));<br>const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs));<br>MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV);<br>return wrap(LT);'<br><br><br>srcsparse <- '<br>using Eigen::Map; <br>using Eigen::SparseMatrix;<br>using Eigen::LLT;<br><br>const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));<br>const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));<br>Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver;<br>SparseMatrix<double> LT = solver.compute(Omega).matrixL()<b>.solve(invSqrtV);</b><br>return wrap(LT);'<br><br></div><div>The problem is that there is no .solve() method associated with the matrixL(). <br></div><div>It is just a intermediary step when computing<br><br>solver.compute(Omega).solve(invSqrtV)<br><br></div><div>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><br></div><div>Thank you.<span class="HOEnZb"><font color="#888888"><br></font></span></div><span class="HOEnZb"><font color="#888888"><div>-- <br><div><div dir="ltr"><div><div dir="ltr"><div>Wagner Hugo Bonat<br>----------------------------------------------------------------------------------------------<br>Department of Mathematics and Computer Science (IMADA)<br>University of Southern Denmark (SDU) and<br>Laboratório de Estatística e Geoinformação (LEG)<br>Universidade Federal do Paraná (UFPR)<br><br></div></div></div></div></div>
</div></font></span></div></div>
<br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">Yixuan Qiu <<a href="mailto:yixuan.qiu@cos.name" target="_blank">yixuan.qiu@cos.name</a>><br>Department of Statistics,<br>Purdue University<br></div>
</div>