qsolve {INLA} | R Documentation |
This routine use the GMRFLib implementation to solve linear systems with a SPD matrix.
inla.qsolve(Q, B, reordering = inla.reorderings(), method = c("solve", "forward", "backward"))
Q |
A SPD matrix, either as a (dense) matrix or sparse-matrix |
.
B |
The right hand side matrix, either as a (dense) matrix or sparse-matrix. |
reordering |
The type of reordering algorithm to be used for |
method |
The system to solve, one of |
inla.qsolve
returns a matrix X
,
which is the solution of Q X = B
, L X = B
or L^T X = B
depending on the value of method
.
Havard Rue hrue@r-inla.org
n = 10 nb <- n-1 QQ = matrix(rnorm(n^2), n, n) QQ <- QQ %*% t(QQ) Q = inla.as.sparse(QQ) B = matrix(rnorm(n*nb), n, nb) X = inla.qsolve(Q, B, method = "solve") XX = inla.qsolve(Q, B, method = "solve", reordering = inla.qreordering(Q)) print(paste("err solve1", sum(abs( Q %*% X - B)))) print(paste("err solve2", sum(abs( Q %*% XX - B)))) ## the forward and backward solve is tricky, as after permutation and with Q=LL', then L is ## lower triangular, but L in the orginal ordering is not lower triangular. if the rhs is iid ## noise, this is not important. to control the reordering, then the 'taucs' library must be ## used. inla.setOption(smtp = 'taucs') ## case 1. use the matrix as is, no reordering r <- "identity" L = t(chol(Q)) X = inla.qsolve(Q, B, method = "forward", reordering = r) XX = inla.qsolve(Q, B, method = "backward", reordering = r) print(paste("err forward ", sum(abs(L %*% X - B)))) print(paste("err backward", sum(abs(t(L) %*% XX - B)))) ## case 2. use a reordering from the library r <- inla.qreordering(Q) im <- r$ireordering m <- r$reordering print(cbind(idx = 1:n, m, im) ) Qr <- Q[im, im] L = t(chol(Qr))[m, m] X = inla.qsolve(Q, B, method = "forward", reordering = r) XX = inla.qsolve(Q, B, method = "backward", reordering = r) print(paste("err forward ", sum(abs( L %*% X - B)))) print(paste("err backward", sum(abs( t(L) %*% XX - B))))