<div dir="ltr"><div><div>Dear all,<br><br></div>As requested I will try provide a better explanation. Consider the code below:<br><div><br></div>Let C = V^{1/2}%*%Omega%*%V^{1/2} and C^{1/2} be the Cholesky decomposition, such that C^{1/2}%*%C^{1/2}^T = C.<br>My goal is to compute Omega^{-1/2}%*%V{-1/2}. <br>My idea is to use that I know V^{-1/2} and rewrite C^{-1} = V^{-1/2}%*%Omega^{-1}%*%V^{-1/2} = V^{-1/2}%*%Omega^{-1/2}%*%(Omega^{-1/2}%*%V^{-1/2})^T<br>The
main term now is Omega^{-1/2}%*%V^{-1/2} for computing it my idea is
first compute the Cholesky denoted by Omega^{1/2} and then using the
triangular structure to solve a system to get Omega^{-1/2}%*%V^{-1/2}.<br><br></div>Some auxiliary functions to provide an example.<br><div><div><br># Inverse of square root of V<br>build.invSqrtV <- function(mu, power,n.sample){Diagonal(n.sample, 1/sqrt(mu^power))}<br><br>## Omega matrix - It is an example of dense structure.<br>build.Omega <- function(par,U){par[1]*exp(-U/par[2])}<br><br>## Small sample size<br>x1 <- seq(0,1,l=30)<br>x2 <- seq(0,1,l=30)<br>grid <- expand.grid(x1,x2)<br>U <- forceSymmetric(as.matrix(dist(grid)))<br>cov <- seq(-1,1,l=900)<br>X <- model.matrix(~ cov)<br>mu <- exp(X%*%c(1,0.5))<br><br>Vinvsqrt <- build.invSqrtV(mu, power = 2, n.sample =900)<br>Omega <- build.Omega(par=c(1,0.25),U)<br><br></div><br><div>## R function<br></div><div>my.partial.solveR <- function(Omega,invSqrtV){ solve(t(chol(Omega)), invSqrtV)}<br>result1 <- my.partial.solveR(Omega, Vinvsqrt)<br><br>## Now using RcppEigen<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>## Here I need to convert from dsyMatrix to matrix<br>Omega.m <- as.matrix(Omega)<br>Vinvsqrt.m <- as.matrix(Vinvsqrt)<br>my.partial.solveC <- cxxfunction(signature(Omegas = "mat", invSqrtVs = "mat"), body = src, plugin = "RcppEigen")<br>result2 <- my.partial.solveC(Omega.m,Vinvsqrt.m)<br><br>## Comparing the results<br>tt <- abs(result1) - abs(result2)<br>sum(tt)<br><br>## Comparing the computational time<br>benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m, Vinvsqrt.m), order = "relative")<br><br></div><div> test replications elapsed relative<br>2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 1.000<br>1 my.partial.solveR(Omega, Vinvsqrt) 100 25.414 2.855<br><br></div><div>You can see that the C++ code is much faster than my R code.<br><br>But, now consider that my Omega matrix is sparse<br>z <- rep(1,10)<br>Z <- z%*%t(z)<br>listZ <- list()<br>for(i in 1:90){listZ[[i]] <- Z}<br>Omega <- Diagonal(900,1) + 0.8*bdiag(listZ)<br>Omega.m <- as.matrix(Omega)<br><br>benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m, Vinvsqrt.m), order = "relative")<br></div><div><br> test replications elapsed relative<br>1 my.partial.solveR(Omega, Vinvsqrt) 100 0.101 1.000<br>2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 88.129<br><div><br></div><div>Now, the results are completely different because my C++ code does not take into account the sparse structure.<br></div><div>My idea is reproduze the same idea using SparseMatrix class. I try the following function, but it does not work, since there is no solve method associated with the solver. So, I do not know how to use the triangular structure to solve the linear system<br><br>src <- '<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::Upper> solver;<br>solver.compute(Omega);<br>SparseMatrix<double> L = solver.matrixL().solve(invSqrtV); #### Do not work<br>return wrap(L);'<br></div><div><br></div><div>Any help is welcome. <br><br></div><div>Thank you!<br></div><div><br><br></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-06-09 7:07 GMT+02:00 Yixuan Qiu <span dir="ltr"><<a href="mailto:yixuan.qiu@cos.name" target="_blank">yixuan.qiu@cos.name</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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"><div><div class="h5">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></div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><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><font color="#888888"><br></font></span></div><span><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></div></div><span class="">_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">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></span></blockquote></div><span class="HOEnZb"><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div>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>
</font></span></div>
</blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature"><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>