RandomFieldsUtils (version 0.5.3)

Cholesky: Cholesky Decomposition of Positive Definite Matrices

Description

This function calculates the Cholesky decomposition of a matrix.

Usage

cholx(a)
chol2mv(C, n)
tcholRHS(C, RHS)

Arguments

a

a square real-valued positive definite matrix

C

a (pivoted) Cholesky decomposition calculated by cholx

n

integer. Number of realisations of the multivariate normal distribution

RHS

vector

Value

cholx returns a matrix containing the Cholesky decomposition (in its upper part).

chol2mv takes the Cholesky decomposition and returns a n realisations of a multivariate normal distribution with mean 0 and covariance function a

tcholRHS multiplies the vector RHS from the right to lower triangular matrix of the Cholesky decomposition. See examples below.

Details

If the matrix is diagonal direct calculations are performed.

Else the Cholesky decomposition is tried.

See Also

chol.spam in the package spam

Examples

# NOT RUN {
if (FALSE) {
## This examples shows that 'cholesky' can be much faster
## than 'chol'

## creating a covariance matrix for a temporal process
covmatrix <- function(model, x) {
  x <- as.matrix(dist(x))
  return(eval(substitute(model)))
}

size <- 600
x <- runif(size, 0, size / 10)
M <- covmatrix((1 - x) * (x < 1) , x) ## Askey's model of covariance
b <- seq(0, 1, len=size)
system.time(C2 <- chol(M))
system.time(C1 <- cholx(M))
range(C2 - C1)
stopifnot(all(abs(C2 - C1) < 10^{-9}))
}



##########################
## Example showing the use of chol2mv and tcholRHS
n <- 10
M <- matrix(nc=n, runif(n^2))
M <- M %*% t(M) + diag(n)
C <- cholx(M)
set.seed(0)
v1 <- chol2mv(C, 1)
set.seed(0)
v2 <- tcholRHS(C, rnorm(n))
stopifnot(all(v1 == v2))



##########################
## The following example shows pivoted Cholesky can be used
## and the pivotation permutation can be transferred to
## subsequent Cholesky decompositions
# }
# NOT RUN {
<!-- %      library(RandomFieldsUtils) -->
# }
# NOT RUN {
set.seed(0)
n <- if (interactive()) 1000 else 100
x <- 1:n
y <- runif(n)
M <- x %*% t(x) + rev(x) %*% t(rev(x)) + y %*% t(y)

## do pivoting
RFoptions(pivot = PIVOT_DO)
print(system.time(C <- cholx(M)))
print(range(crossprod(C) - M))
str(C)

## use the same pivoted decomposition as in the previous decomposition
M2 <- M +  n * diag(1:n)
RFoptions(pivot = PIVOT_IDX,
          pivot_idx = attr(C, "pivot_idx"),
          pivot_actual_size = attr(C, "pivot_actual_size"))
print(system.time(C2 <- cholx(M2)))
print(range(crossprod(C2) - M2))
range((crossprod(C2) - M2) / M2)
str(C2)

# }
# NOT RUN {
# }