Sparse Cholesky Factorizations (original) (raw)

Description

CHMfactor is the virtual class of sparse Cholesky factorizations of n \times n real, symmetric matrices A, having the general form

P_1 A P_1' = L_1 D L_1' \overset{D_{jj} \ge 0}{=} L L'

or (equivalently)

A = P_1' L_1 D L_1' P_1 \overset{D_{jj} \ge 0}{=} P_1' L L' P_1

whereP_1 is a permutation matrix,L_1 is a unit lower triangular matrix,D is a diagonal matrix, andL = L_1 \sqrt{D}. The second equalities hold only for positive semidefinite A, for which the diagonal entries of D are non-negative and \sqrt{D} is well-defined.

The implementation of class CHMfactor is based on CHOLMOD's C-level cholmod_factor_struct. Virtual subclasses CHMsimpl and CHMsuper separate the simplicial and supernodal variants. These have nonvirtual subclasses [dn]CHMsimpl and [dn]CHMsuper, where prefix ‘⁠d⁠’ and prefix ‘⁠n⁠’ are reserved for numeric and symbolic factorizations, respectively.

Usage

isLDL(x)

Arguments

Value

isLDL(x) returns TRUE or FALSE:TRUE if x stores the lower triangular entries of L_1-I+D,FALSE if x stores the lower triangular entries of L.

Slots

Of CHMfactor:

Dim, Dimnames

inherited from virtual classMatrixFactorization.

colcount

an integer vector of length Dim[1]giving an estimate of the number of nonzero entries in each column of the lower triangular Cholesky factor. If symbolic analysis was performed prior to factorization, then the estimate is exact.

perm

a 0-based integer vector of length Dim[1]specifying the permutation applied to the rows and columns of the factorized matrix. perm of length 0 is valid and equivalent to the identity permutation, implying no pivoting.

type

an integer vector of length 6 specifying details of the factorization. The elements correspond to members ordering, is_ll, is_super,is_monotonic, maxcsize, and maxesizeof the original cholmod_factor_struct. Simplicial and supernodal factorizations are distinguished by is_super. Simplicial factorizations do not usemaxcsize or maxesize. Supernodal factorizations do not use is_ll or is_monotonic.

Of CHMsimpl (all unused by nCHMsimpl):

nz

an integer vector of length Dim[1]giving the number of nonzero entries in each column of the lower triangular Cholesky factor. There is at least one nonzero entry in each column, because the diagonal elements of the factor are stored explicitly.

p

an integer vector of length Dim[1]+1. Row indices of nonzero entries in column j of the lower triangular Cholesky factor are obtained asi[p[j]+seq_len(nz[j])]+1.

i

an integer vector of length greater than or equal to sum(nz) containing the row indices of nonzero entries in the lower triangular Cholesky factor. These are grouped by column and sorted within columns, but the columns themselves need not be ordered monotonically. Columns may be overallocated, i.e., the number of elements of i reserved for columnj may exceed nz[j].

prv, nxt

integer vectors of lengthDim[1]+2 indicating the order in which the columns of the lower triangular Cholesky factor are stored in iand x. Starting from j <- Dim[1]+2, the recursion j <- nxt[j+1]+1 traverses the columns in forward order and terminates when nxt[j+1] = -1. Starting from j <- Dim[1]+1, the recursion j <- prv[j+1]+1 traverses the columns in backward order and terminates when prv[j+1] = -1.

Of dCHMsimpl:

x

a numeric vector parallel to i containing the corresponding nonzero entries of the lower triangular Cholesky factor L or (if and only if type[2]is 0) of the lower triangular matrix L_1-I+D.

Of CHMsuper:

super, pi, px

integer vectors of length nsuper+1, where nsuper is the number of supernodes. super[j]+1 is the index of the leftmost column of supernode j. The row indices of supernodej are obtained as s[pi[j]+seq_len(pi[j+1]-pi[j])]+1. The numeric entries of supernode j are obtained asx[px[j]+seq_len(px[j+1]-px[j])]+1 (if slot xis available).

s

an integer vector of length greater than or equal to Dim[1] containing the row indices of the supernodes.s may contain duplicates, but not within a supernode, where the row indices must be increasing.

Of dCHMsuper:

x

a numeric vector of length less than or equal toprod(Dim) containing the numeric entries of the supernodes.

Extends

Class MatrixFactorization, directly.

Instantiation

Objects can be generated directly by calls of the formnew("dCHMsimpl", ...), etc., but dCHMsimpl anddCHMsuper are more typically obtained as the value ofCholesky(x, ...) for x inheriting fromsparseMatrix(often dsCMatrix).

There is currently no API outside of calls to newfor generating nCHMsimpl and nCHMsuper. These classes are vestigial and may be formally deprecated in a future version of Matrix.

Methods

coerce

signature(from = "CHMsimpl", to = "dtCMatrix"): returns a dtCMatrix representing the lower triangular Cholesky factor L _or_the lower triangular matrix L_1-I+D, the latter if and only if from@type[2] is 0.

coerce

signature(from = "CHMsuper", to = "dgCMatrix"): returns a dgCMatrix representing the lower triangular Cholesky factor L. Note that, for supernodes spanning two or more columns, the supernodal algorithm by design stores non-structural zeros above the main diagonal, hence dgCMatrix is indeed more appropriate than dtCMatrixas a coercion target.

determinant

signature(from = "CHMfactor", logarithm = "logical"): behaves according to an optional argument sqrt. If sqrt = FALSE, then this method computes the determinant of the factorized matrix A or its logarithm. If sqrt = TRUE, then this method computes the determinant of the factor L = L_1 sqrt(D) or its logarithm, giving NaN for the modulus when Dhas negative diagonal elements. For backwards compatibility, the default value of sqrt is TRUE, but that can be expected change in a future version of Matrix, hence defensive code will always set sqrt (to TRUE, if the code must remain backwards compatible with Matrix < 1.6-0). Calls to this method not setting sqrtmay warn about the pending change. The warnings can be disabled with options(Matrix.warnSqrtDefault = 0).

diag

signature(x = "CHMfactor"): returns a numeric vector of length n containing the diagonal elements of D, which (if they are all non-negative) are the squared diagonal elements of L.

expand

signature(x = "CHMfactor"): see expand-methods.

expand1

signature(x = "CHMsimpl"): see expand1-methods.

expand1

signature(x = "CHMsuper"): see expand1-methods.

expand2

signature(x = "CHMsimpl"): see expand2-methods.

expand2

signature(x = "CHMsuper"): see expand2-methods.

image

signature(x = "CHMfactor"): see image-methods.

nnzero

signature(x = "CHMfactor"): see nnzero-methods.

solve

signature(a = "CHMfactor", b = .): see solve-methods.

update

signature(object = "CHMfactor"): returns a copy of object with the same nonzero pattern but with numeric entries updated according to additional arguments parent and mult, where parentis (coercible to) a dsCMatrix or adgCMatrix and mult is a numeric vector of positive length.
The numeric entries are updated with those of the Cholesky factor of F(parent) + mult[1] * I, i.e.,F(parent) plus mult[1] times the identity matrix, where F = identity for symmetric parentand F = tcrossprod for other parent. The nonzero pattern of F(parent) must match that of S if object = Cholesky(S, ...).

updown

signature(update = ., C = ., object = "CHMfactor"): see updown-methods.

References

The CHOLMOD source code; seehttps://github.com/DrTimothyAldenDavis/SuiteSparse, notably the header file ‘CHOLMOD/Include/cholmod.h’ defining cholmod_factor_struct.

Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate.ACM Transactions on Mathematical Software,35(3), Article 22, 1-14. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1145/1391989.1391995")}

Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004). Algorithm 837: AMD, an approximate minimum degree ordering algorithm.ACM Transactions on Mathematical Software,17(4), 886-905. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1145/1024074.1024081")}

Golub, G. H., & Van Loan, C. F. (2013).Matrix computations (4th ed.). Johns Hopkins University Press. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.56021/9781421407944")}

See Also

Class dsCMatrix.

Generic functions Cholesky, updown,expand1 and expand2.

Examples

showClass("dCHMsimpl")
showClass("dCHMsuper")
set.seed(2)

m <- 1000L
n <- 200L
M <- rsparsematrix(m, n, 0.01)
A <- crossprod(M)

## With dimnames, to see that they are propagated :
dimnames(A) <- dn <- rep.int(list(paste0("x", seq_len(n))), 2L)

(ch.A <- Cholesky(A)) # pivoted, by default
str(e.ch.A <- expand2(ch.A, LDL =  TRUE), max.level = 2L)
str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L)

ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...)
ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...)

## A ~ P1' L1 D L1' P1 ~ P1' L L' P1 in floating point
stopifnot(exprs = {
    identical(names(e.ch.A), c("P1.", "L1", "D", "L1.", "P1"))
    identical(names(E.ch.A), c("P1.", "L" ,      "L." , "P1"))
    identical(e.ch.A[["P1"]],
              new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
                  margin = 2L, perm = invertPerm(ch.A@perm, 0L, 1L)))
    identical(e.ch.A[["P1."]], t(e.ch.A[["P1"]]))
    identical(e.ch.A[["L1."]], t(e.ch.A[["L1"]]))
    identical(E.ch.A[["L." ]], t(E.ch.A[["L" ]]))
    identical(e.ch.A[["D"]], Diagonal(x = diag(ch.A)))
    all.equal(E.ch.A[["L"]], with(e.ch.A, L1 %*% sqrt(D)))
    ae1(A, with(e.ch.A, P1. %*% L1 %*% D %*% L1. %*% P1))
    ae1(A, with(E.ch.A, P1. %*% L  %*%         L.  %*% P1))
    ae2(A[ch.A@perm + 1L, ch.A@perm + 1L], with(e.ch.A, L1 %*% D %*% L1.))
    ae2(A[ch.A@perm + 1L, ch.A@perm + 1L], with(E.ch.A, L  %*%         L. ))
})

## Factorization handled as factorized matrix
## (in some cases only optionally, depending on arguments)
b <- rnorm(n)
stopifnot(identical(det(A), det(ch.A, sqrt = FALSE)),
          identical(solve(A, b), solve(ch.A, b, system = "A")))

u1 <- update(ch.A,   A , mult = sqrt(2))
u2 <- update(ch.A, t(M), mult = sqrt(2)) # updating with crossprod(M), not M
stopifnot(all.equal(u1, u2, tolerance = 1e-14))