Dear all,
As requested I will try provide a better explanation. Consider the code
below:
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.
My goal is to compute Omega^{-1/2}%*%V{-1/2}.
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
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}.
Some auxiliary functions to provide an example.
# Inverse of square root of V
build.invSqrtV <- function(mu, power,n.sample){Diagonal(n.sample,
1/sqrt(mu^power))}
## Omega matrix - It is an example of dense structure.
build.Omega <- function(par,U){par[1]*exp(-U/par[2])}
## Small sample size
x1 <- seq(0,1,l=30)
x2 <- seq(0,1,l=30)
grid <- expand.grid(x1,x2)
U <- forceSymmetric(as.matrix(dist(grid)))
cov <- seq(-1,1,l=900)
X <- model.matrix(~ cov)
mu <- exp(X%*%c(1,0.5))
Vinvsqrt <- build.invSqrtV(mu, power = 2, n.sample =900)
Omega <- build.Omega(par=c(1,0.25),U)
## R function
my.partial.solveR <- function(Omega,invSqrtV){ solve(t(chol(Omega)),
invSqrtV)}
result1 <- my.partial.solveR(Omega, Vinvsqrt)
## Now using RcppEigen
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);'
## Here I need to convert from dsyMatrix to matrix
Omega.m <- as.matrix(Omega)
Vinvsqrt.m <- as.matrix(Vinvsqrt)
my.partial.solveC <- cxxfunction(signature(Omegas = "mat", invSqrtVs =
"mat"), body = src, plugin = "RcppEigen")
result2 <- my.partial.solveC(Omega.m,Vinvsqrt.m)
## Comparing the results
tt <- abs(result1) - abs(result2)
sum(tt)
## Comparing the computational time
benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m,
Vinvsqrt.m), order = "relative")
test replications elapsed relative
2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 1.000
1 my.partial.solveR(Omega, Vinvsqrt) 100 25.414 2.855
You can see that the C++ code is much faster than my R code.
But, now consider that my Omega matrix is sparse
z <- rep(1,10)
Z <- z%*%t(z)
listZ <- list()
for(i in 1:90){listZ[[i]] <- Z}
Omega <- Diagonal(900,1) + 0.8*bdiag(listZ)
Omega.m <- as.matrix(Omega)
benchmark(my.partial.solveR(Omega, Vinvsqrt), my.partial.solveC(Omega.m,
Vinvsqrt.m), order = "relative")
test replications elapsed relative
1 my.partial.solveR(Omega, Vinvsqrt) 100 0.101 1.000
2 my.partial.solveC(Omega.m, Vinvsqrt.m) 100 8.901 88.129
Now, the results are completely different because my C++ code does not take
into account the sparse structure.
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
src <- '
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::Upper> solver;
solver.compute(Omega);
SparseMatrix<double> L = solver.matrixL().solve(invSqrtV); #### Do not work
return wrap(L);'
Any help is welcome.
Thank you!
2015-06-09 7:07 GMT+02:00 Yixuan Qiu <yixuan....@cos.name>:
> 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
>
--
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
