CHMfactor-class: Sparse Cholesky Factorizations
In Matrix: Sparse and Dense Matrix Classes and Methods
CHMfactor-class R Documentation
Sparse Cholesky Factorizations
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
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, Dimnamesinherited from virtual class
MatrixFactorization.
colcountan 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.
perma 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.
typean 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):
nzan 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.
pan 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.
ian 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, nxtinteger 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:
xa 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, pxinteger 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).
san 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:
xa 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
coercesignature(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.
coercesignature(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.
determinantsignature(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).
diagsignature(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.
expandsignature(x = "CHMfactor"):
see expand-methods.
expand1signature(x = "CHMsimpl"):
see expand1-methods.
expand1signature(x = "CHMsuper"):
see expand1-methods.
expand2signature(x = "CHMsimpl"):
see expand2-methods.
expand2signature(x = "CHMsuper"):
see expand2-methods.
imagesignature(x = "CHMfactor"):
see image-methods.
nnzerosignature(x = "CHMfactor"):
see nnzero-methods.
solvesignature(a = "CHMfactor", b = .):
see solve-methods.
updatesignature(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, ...).
updownsignature(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.
\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))
Matrix documentation built on Oct. 19, 2024, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| CHMfactor-class | R Documentation |
Sparse Cholesky Factorizations
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
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 |
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,Dimnamesinherited from virtual class
MatrixFactorization.colcountan 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.perma 0-based integer vector of length
Dim[1]specifying the permutation applied to the rows and columns of the factorized matrix.permof length 0 is valid and equivalent to the identity permutation, implying no pivoting.typean integer vector of length 6 specifying details of the factorization. The elements correspond to members
ordering,is_ll,is_super,is_monotonic,maxcsize, andmaxesizeof the originalcholmod_factor_struct. Simplicial and supernodal factorizations are distinguished byis_super. Simplicial factorizations do not usemaxcsizeormaxesize. Supernodal factorizations do not useis_lloris_monotonic.
Of CHMsimpl (all unused by nCHMsimpl):
nzan 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.pan integer vector of length
Dim[1]+1. Row indices of nonzero entries in columnjof the lower triangular Cholesky factor are obtained asi[p[j]+seq_len(nz[j])]+1.ian 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 ofireserved for columnjmay exceednz[j].prv,nxtinteger vectors of length
Dim[1]+2indicating the order in which the columns of the lower triangular Cholesky factor are stored iniandx. Starting fromj <- Dim[1]+2, the recursionj <- nxt[j+1]+1traverses the columns in forward order and terminates whennxt[j+1] = -1. Starting fromj <- Dim[1]+1, the recursionj <- prv[j+1]+1traverses the columns in backward order and terminates whenprv[j+1] = -1.
Of dCHMsimpl:
xa numeric vector parallel to
icontaining the corresponding nonzero entries of the lower triangular Cholesky factorLor (if and only iftype[2]is 0) of the lower triangular matrixL_1-I+D.
Of CHMsuper:
super,pi,pxinteger vectors of length
nsuper+1, wherensuperis the number of supernodes.super[j]+1is the index of the leftmost column of supernodej. The row indices of supernodejare obtained ass[pi[j]+seq_len(pi[j+1]-pi[j])]+1. The numeric entries of supernodejare obtained asx[px[j]+seq_len(px[j+1]-px[j])]+1(if slotxis available).san integer vector of length greater than or equal to
Dim[1]containing the row indices of the supernodes.smay contain duplicates, but not within a supernode, where the row indices must be increasing.
Of dCHMsuper:
xa 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
coercesignature(from = "CHMsimpl", to = "dtCMatrix"): returns adtCMatrixrepresenting the lower triangular Cholesky factorLor the lower triangular matrixL_1-I+D, the latter if and only iffrom@type[2]is 0.coercesignature(from = "CHMsuper", to = "dgCMatrix"): returns adgCMatrixrepresenting the lower triangular Cholesky factorL. Note that, for supernodes spanning two or more columns, the supernodal algorithm by design stores non-structural zeros above the main diagonal, hencedgCMatrixis indeed more appropriate thandtCMatrixas a coercion target.determinantsignature(from = "CHMfactor", logarithm = "logical"): behaves according to an optional argumentsqrt. Ifsqrt = FALSE, then this method computes the determinant of the factorized matrixAor its logarithm. Ifsqrt = TRUE, then this method computes the determinant of the factorL = L_1 sqrt(D)or its logarithm, givingNaNfor the modulus whenDhas negative diagonal elements. For backwards compatibility, the default value ofsqrtisTRUE, but that can be expected change in a future version of Matrix, hence defensive code will always setsqrt(toTRUE, if the code must remain backwards compatible with Matrix< 1.6-0). Calls to this method not settingsqrtmay warn about the pending change. The warnings can be disabled withoptions(Matrix.warnSqrtDefault = 0).diagsignature(x = "CHMfactor"): returns a numeric vector of lengthncontaining the diagonal elements ofD, which (if they are all non-negative) are the squared diagonal elements ofL.expandsignature(x = "CHMfactor"): seeexpand-methods.expand1signature(x = "CHMsimpl"): seeexpand1-methods.expand1signature(x = "CHMsuper"): seeexpand1-methods.expand2signature(x = "CHMsimpl"): seeexpand2-methods.expand2signature(x = "CHMsuper"): seeexpand2-methods.imagesignature(x = "CHMfactor"): seeimage-methods.nnzerosignature(x = "CHMfactor"): seennzero-methods.solvesignature(a = "CHMfactor", b = .): seesolve-methods.updatesignature(object = "CHMfactor"): returns a copy ofobjectwith the same nonzero pattern but with numeric entries updated according to additional argumentsparentandmult, whereparentis (coercible to) adsCMatrixor adgCMatrixandmultis a numeric vector of positive length.
The numeric entries are updated with those of the Cholesky factor ofF(parent) + mult[1] * I, i.e.,F(parent)plusmult[1]times the identity matrix, whereF = identityfor symmetricparentandF = tcrossprodfor otherparent. The nonzero pattern ofF(parent)must match that ofSifobject = Cholesky(S, ...).updownsignature(update = ., C = ., object = "CHMfactor"): seeupdown-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. \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))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.