CHMfactor-class function

Sparse Cholesky Factorizations

CHMfactor is the virtual class of sparse Cholesky factorizations of $n-by-n$ real, symmetric matrices $A$, having the general form [REMOVE_ME]$P_1 A P_1' = L_1 D L_1' \overset{D_{jj} \ge 0}{=} L L'P1 * A * P1' = L1 * D * L1' [ = L * L' ] [REMOVE_ME_2]$

or (equivalently) [REMOVE_ME]$A = P_1' L_1 D L_1' P_1 \overset{D_{jj} \ge 0}{=} P_1' L L' P_1A = P1' L1 * D * L1' * P1 [ = P1' * L * L' * P1 ] [REMOVE_ME_2]$

where $P1$ is a permutation matrix, $L1$ is a unit lower triangular matrix, $D$ is a diagonal matrix, and $L = L1 * 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.

class

Description

CHMfactor is the virtual class of sparse Cholesky factorizations of $n-by-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'P1 * A * P1' = L1 * D * L1' [ = L * L' ]$$

or (equivalently)

$$A = P_1' L_1 D L_1' P_1 \overset{D_{jj} \ge 0}{=} P_1' L L' P_1A = P1' L1 * D * L1' * P1 [ = P1' * L * L' * P1 ]$$

where $P1$ is a permutation matrix, $L1$ is a unit lower triangular matrix, $D$ is a diagonal matrix, and $L = L1 * 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.

isLDL(x)

Arguments

• x: an object inheriting from virtual class CHMfactor, almost always the result of a call to generic function Cholesky.

Returns

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

Slots

Of CHMfactor:

• Dim, Dimnames: inherited from virtual class MatrixFactorization.

• 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 maxesize

   of the original `cholmod_factor_struct`. Simplicial and supernodal factorizations are distinguished by `is_super`. Simplicial factorizations do not use `maxcsize` 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 as i[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 column j may exceed nz[j].

• prv, nxt: integer vectors of length Dim[1]+2 indicating the order in which the columns of the lower triangular Cholesky factor are stored in i

   and `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 $L1-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 supernode j are obtained as s[pi[j]+seq_len(pi[j+1]-pi[j])]+1. The numeric entries of supernode j are obtained as x[px[j]+seq_len(px[j+1]-px[j])]+1 (if slot x

   is 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 to prod(Dim) containing the numeric entries of the supernodes.

Extends

Class MatrixFactorization, directly.

Instantiation

Objects can be generated directly by calls of the form new("dCHMsimpl", ...), etc., but dCHMsimpl and dCHMsuper are more typically obtained as the value of Cholesky(x, ...) for x inheriting from sparseMatrix

(often dsCMatrix).

There is currently no API outside of calls to new

for 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 $L1-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 dtCMatrix

   as 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 = L1 * sqrt(D)$ or its logarithm, giving NaN for the modulus when $D$

   has 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 `sqrt`
   
   may 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 parent

   is (coercible to) a `dsCMatrix` or a `dgCMatrix` 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 `parent`
   
   and `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.

See Also

Class dsCMatrix.

Generic functions Cholesky, updown, expand1 and expand2.

References

The CHOLMOD source code; see https://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. 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. tools:::Rd_expr_doi("10.1145/1024074.1024081")

Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. tools:::Rd_expr_doi("10.56021/9781421407944")

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