cholPosDef: Cholesky Decomposition of Positive Definite Matrices

In RandomFieldsUtils: Utilities for the Simulation and Analysis of Random Fields and Genetic Data

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

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.

References

Harbrecht, H., Peters, M., Schneider, R. (2012) On the low-rank approximation by the pivoted Cholesky decomposition. Appl. Num. Math. 62, 428–440.

Examples

##########################
## 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

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, la_mode=LA_INTERN)
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, la_mode=LA_INTERN,
          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)

RFoptions(pivot = PIVOT_AUTO, la_mode = LA_AUTO)

RandomFieldsUtils documentation built on April 19, 2022, 5:09 p.m.