Sparse Cholesky Factorizations

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$$ where $P_1$ is a permutation matrix, $L_1$ is a unit lower triangular matrix, $D$ is a diagonal matrix, and $L = 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

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

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 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 $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 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 $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 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 = L_1 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.

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. 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. doi:10.1145/1024074.1024081

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

Examples

showClass("dCHMsimpl")
#> Class "dCHMsimpl" [package "Matrix"]
#> 
#> Slots:
#>                                                                               
#> Name:         x        p        i       nz      nxt      prv     type colcount
#> Class:  numeric  integer  integer  integer  integer  integer  integer  integer
#>                                  
#> Name:      perm      Dim Dimnames
#> Class:  integer  integer     list
#> 
#> Extends: 
#> Class "CHMsimpl", directly
#> Class "CHMfactor", by class "CHMsimpl", distance 2
#> Class "CholeskyFactorization", by class "CHMsimpl", distance 3
#> Class "MatrixFactorization", by class "CHMsimpl", distance 4
showClass("dCHMsuper")
#> Class "dCHMsuper" [package "Matrix"]
#> 
#> Slots:
#>                                                                               
#> Name:         x    super       pi       px        s     type colcount     perm
#> Class:  numeric  integer  integer  integer  integer  integer  integer  integer
#>                         
#> Name:       Dim Dimnames
#> Class:  integer     list
#> 
#> Extends: 
#> Class "CHMsuper", directly
#> Class "CHMfactor", by class "CHMsuper", distance 2
#> Class "CholeskyFactorization", by class "CHMsuper", distance 3
#> Class "MatrixFactorization", by class "CHMsuper", distance 4
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
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#>   ..@ x       : num [1:10105] 4.2231 0.0947 0.0831 -0.2842 -0.1314 ...
#>   ..@ p       : int [1:201] 0 10 20 43 54 66 76 84 105 138 ...
#>   ..@ i       : int [1:10105] 0 102 107 112 122 159 161 177 191 195 ...
#>   ..@ nz      : int [1:200] 10 10 23 11 12 10 8 21 33 6 ...
#>   ..@ nxt     : int [1:202] 1 2 3 4 5 6 7 8 9 10 ...
#>   ..@ prv     : int [1:202] 201 0 1 2 3 4 5 6 7 8 ...
#>   ..@ type    : int [1:6] 2 0 0 1 0 0
#>   ..@ colcount: int [1:200] 10 10 23 11 12 10 8 21 33 6 ...
#>   ..@ perm    : int [1:200] 106 60 28 149 188 84 178 102 98 185 ...
#>   ..@ Dim     : int [1:2] 200 200
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:200] "x1" "x2" "x3" "x4" ...
#>   .. ..$ : chr [1:200] "x1" "x2" "x3" "x4" ...
str(e.ch.A <- expand2(ch.A, LDL =  TRUE), max.level = 2L)
#> List of 5
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#>  $ L1 :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ D  :Formal class 'ddiMatrix' [package "Matrix"] with 4 slots
#>  $ L1.:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ P1 :Formal class 'pMatrix' [package "Matrix"] with 4 slots
str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L)
#> List of 4
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#>  $ L  :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ L. :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ P1 :Formal class 'pMatrix' [package "Matrix"] with 4 slots

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