spam (version 0.13-3)

chol: Cholesky Factorization for Sparse Matrices

Description

chol performs a Cholesky decomposition of a symmetric positive definite sparse matrix x of class spam.

Usage

chol(x, ...)

Arguments

x

symmetric positive definite matrix of class spam.

...

further arguments passed to or from other methods.

Value

• The function returns the Cholesky factor in an object of class spam.chol.method. Recall that the latter is the Cholesky factor of a reordered matrix x, see also ordering.

Details

chol performs a Cholesky decomposition of a symmetric positive definite sparse matrix x of class spam. Currently, there is only the block sparse Cholesky algorithm of Ng and Peyton (1993) implemented (method=NgPeyton), based on a minimum degree permutation (ordering=MinDeg). The Cholesky routine is essentially identitcal to what is available in the package SparseM. The output of chol can be used with forwardsolve backsolve to solve a system of linear equations. The Cholesky decompositions may require additional parameters, linked to memory allocation. These memory parameters and their default values for the NgPeyton method are tmpmax=500*nrow(x), nsubmax=length(x), nnzlmax=max(4*nsubmax, floor(.2*nsubmax^1.3)). In practice, only the first one needs to be set to a larger value if one calculates the Cholesky factor from a virtually full matrix.

References

Ng, E. G. and B. W. Peyton (1993), "Block sparse Cholesky algorithms on advanced uniprocessor computers", SIAM J. Sci. Comput., 14, pp. 1034-1056.

Koenker, R and Ng, P. (2002). SparseM: A Sparse Matrix Package for R, http://www.econ.uiuc.edu/~roger/research

See Also

det, solve, forwardsolve, backsolve and ordering.

Examples

# generate multivariate normals:
n <- 25    # dimension
N <- 1000  # sample size
Sigma <- .25^abs(outer(1:n,1:n,"-"))
Sigma <- as.spam( Sigma, eps=1e-4)

cholS <- chol( Sigma)    
# cholS is the upper triangular part of the permutated matrix Sigma
iord <- ordering(cholS, inv=TRUE)

R <- as.spam(cholS)[,iord]
set.seed(13)
mvsample <- array(rnorm(N*n),c(N,n)) %*% R

# 'mvsample' is of class 'spam'. We need to transform it to a
# regular matrix, as there is no method 'var' for 'spam' (should there?).
norm( var( as.matrix( mvsample)) - Sigma, type="HS")
norm( t(R) %*% R - Sigma, type="sup")