Cholesky

Cholesky Decomposition of Positive Definite Matrices

This function calculates the Cholesky decomposition of a matrix.

Keywords
math
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

Details

If the matrix is diagonal direct calculations are performed.

Else the Cholesky decomposition is tried.

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.

See Also

chol.spam in the package spam

Aliases
  • cholesky
  • chol
  • cholx
  • cholPosDef
  • chol2mv
  • tcholRHS
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 {
# }
Documentation reproduced from package RandomFieldsUtils, version 0.5.3, License: GPL (>= 3)

Community examples

Looks like there are no examples yet.