chol performs a Cholesky
decomposition of a symmetric positive definite sparse matrix
1 2 3 4 5 6 7 8 9
symmetric positive definite matrix of class
should the matrix be permuted, and if, with what algorithm, see ‘Details’ below.
Parameters specific to the method, see ‘Details’ below.
threshold to test symmetry. Defaults to
further arguments passed to or from other methods.
an object from a previous call to
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 (
To pivot/permute the matrix, you can choose between the multiple minimum
pivot="MMD") or reverse Cuthill-Mckee (
from George and Lui (1981). It is also possible to furnish a specific
permutation in which case
pivot is a vector. For compatibility
pivot can also take a logical in which for
no permutation is done and for
TRUE is equivalent to
Often the sparsity structure is fixed and does not change, but the
entries do. In those cases, we can update the Cholesky factor with
update.spam.chol.NgPeyton by suppling a Cholesky factor and the
updated matrix. Notice that the structure is effectively
update(object, x). The update feature without assignement has been disabled.
cholupdatesingular determines how singular matrices
are handled by
update. The function hands back an error
"error"), a warning (
"warning") or the value
The Cholesky decompositions requires parameters, linked to memory
allocation. If the default values are too small the Fortran routine
returns an error to R, which allocates more space and calls the Fortran
routine again. The user can also pass better estimates of the allocation
chol with the argument
nnzcolindices=...). The minimal sizes for a fixed sparsity
structure can be obtained from a
summary call, see ‘Examples’.
The output of
chol can be used with
backsolve to solve a system of linear equations.
Notice that the Cholesky factorization of the package
SparseM is also
based on the algorithm of Ng and Peyton (1993). Whereas the Cholesky
routine of the package
Matrix are based on
CHOLMOD by Timothy A. Davis (
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
Although the symmetric structure of
x is needed, only the upper
diagonal entries are used. By default, the code does check for
symmetry (contrarily to
depending on the matrix size, this is a time consuming test.
A test is ignored if
spam.options( "cholsymmetrycheck") is set to
If a permutation is supplied with
spam.options( "cholpivotcheck") determines if the permutation is
tested for validity (defaults to
Reinhard Furrer, based on Ng and Peyton (1993) Fortran routines
Ng, E. G. and Peyton, B. W. (1993) Block sparse Cholesky algorithms on advanced uniprocessor computers, SIAM J. Sci. Comput., 14, 1034–1056.
Gilbert, J. R., Ng, E. G. and Peyton, B. W. (1994) An efficient algorithm to compute row and column counts for sparse Cholesky factorization, SIAM J. Matrix Anal. Appl., 15, 1075–1091.
George, A. and Liu, J. (1981) Computer Solution of Large Sparse Positive Definite Systems, Prentice Hall.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
# generate multivariate normals: set.seed(13) 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) mvsample <- ( array(rnorm(N*n),c(N,n)) %*% R)[,iord] # It is often better to order the sample than the matrix # R itself. # '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='m') norm( t(R) %*% R - Sigma) # To speed up factorizations, memory allocations can be optimized: opt <- summary(cholS) # here, some elements of Sigma may be changed... cholS <- chol( Sigma, memory=list(nnzR=opt$nnzR,nnzcolindices=opt$nnzc))