Solves linear SPD systems

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")
)

Arguments

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 TAUCS; either one of the names listed in inla.reorderings() or the output from inla.qreordering(Q). The default is "auto" which try several reordering algorithm and use the best one for this particular matrix (using the TAUCS library).

method

The system to solve, one of "solve", "forward" or "backward". Let Q = L L^T, where L is lower triangular (the Cholesky triangle), then method="solve" solves L L^T X = B or equivalently Q X = B, method="forward" solves L X = B, and method="backward" solves L^T X = B.

Value

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.

Author

Havard Rue hrue@r-inla.org

Examples

 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))))
#> [1] "err solve1 6.88615831023753e-14"
 print(paste("err solve2", sum(abs( Q %*% XX - B))))
#> [1] "err solve2 6.88615831023753e-14"

 ## 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))))
#> [1] "err forward  2.60208521396521e-14"
 print(paste("err backward", sum(abs(t(L) %*% XX - B))))
#> [1] "err backward 1.8846035843012e-14"

 ## 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) )
#>       idx  m im
#>  [1,]   1  1  1
#>  [2,]   2  9  9
#>  [3,]   3  8  8
#>  [4,]   4  7  7
#>  [5,]   5  6  6
#>  [6,]   6  5  5
#>  [7,]   7  4  4
#>  [8,]   8  3  3
#>  [9,]   9  2  2
#> [10,]  10 10 10
 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))))
#> [1] "err forward  2.49106291150269e-14"
 print(paste("err backward", sum(abs( t(L) %*% XX - B))))
#> [1] "err backward 3.00592883917261e-14"