The virtual class "CHMfactor" is a class of
CHOLMOD-based Cholesky factorizations of symmetric, sparse,
compressed, column-oriented matrices. Such a factorization is
simplicial (virtual class "CHMsimpl") or supernodal (virtual
class "CHMsuper"). Objects that inherit from these classes are
either numeric factorizations (classes "dCHMsimpl" and
"dCHMsuper") or symbolic factorizations (classes
"nCHMsimpl" and "nCHMsuper").
isLDL(x)# S4 method for CHMfactor
update(object, parent, mult = 0, ...)
.updateCHMfactor(object, parent, mult)
## and many more methods, notably,
## solve(a, b, system = c("A","LDLt","LD","DLt","L","Lt","D","P","Pt"), ...)
## ----- see below
a "CHMfactor" object (almost always the result
of Cholesky()).
a "dsCMatrix" or
"dgCMatrix" matrix object with the same nonzero
pattern as the matrix that generated object. If parent
is symmetric, of class "dsCMatrix", then
object should be a decomposition of a matrix with the same
nonzero pattern as parent. If parent is not symmetric then object
should be the decomposition of a matrix with the same nonzero
pattern as tcrossprod(parent).
Since Matrix version 1.0-8, other "sparseMatrix" matrices are
coerced to dsparseMatrix and
CsparseMatrix if needed.
a numeric scalar (default 0). mult times the
identity matrix is (implicitly) added to parent or
tcrossprod(parent) before updating the decomposition
object.
potentially further arguments to the methods.
Objects can be created by calls of the form new("dCHMsuper", ...)
but are more commonly created via Cholesky(),
applied to dsCMatrix or
lsCMatrix objects.
For an introduction, it may be helpful to look at the expand()
method and examples below.
of "CHMfactor" and all classes inheriting from it:
perm:An integer vector giving the 0-based permutation of the rows and columns chosen to reduce fill-in and for post-ordering.
colcount:Object of class "integer" ....
type:Object of class "integer" ....
Slots of the non virtual classes “[dl]CHM(super|simpl)”:
p:Object of class "integer" of pointers, one
for each column, to the initial (zero-based) index of elements in
the column. Only present in classes that contain "CHMsimpl".
i:Object of class "integer" of length nnzero
(number of non-zero elements). These are the row numbers for
each non-zero element in the matrix. Only present in classes that
contain "CHMsimpl".
x:For the "d*" classes: "numeric" - the
non-zero elements of the matrix.
(x) returns a logical indicating if
x is an \(LDL'\) decomposition or (when FALSE) an
\(LL'\) one.
signature(from = "CHMfactor", to = "sparseMatrix")
(or equivalently, to = "Matrix" or to = "triangularMatrix")
as(*, "sparseMatrix") returns the lower triangular factor
\(L\) from the \(LL'\) form of the Cholesky factorization.
Note that (currently) the factor from the \(LL'\) form is always
returned, even if the "CHMfactor" object represents an
\(LDL'\) decomposition.
Furthermore, this is the factor after any fill-reducing
permutation has been applied. See the expand method for
obtaining both the permutation matrix, \(P\), and the lower
Cholesky factor, \(L\).
signature(from = "CHMfactor", to = "pMatrix")
returns the permutation matrix \(P\), representing the
fill-reducing permutation used in the decomposition.
signature(x = "CHMfactor") returns a list with
components P, the matrix representing the fill-reducing
permutation, and L, the lower triangular Cholesky factor.
The original positive-definite matrix \(A\) corresponds to the product
\(A = P'LL'P\). Because of fill-in during the decomposition the
product may apparently have more non-zeros than the original
matrix, even after applying drop0 to it. However,
the extra "non-zeros" should be very small in magnitude.
signature(x = "CHMfactor"): Plot the image of the
lower triangular factor, \(L\), from the decomposition. This method
is equivalent to image(as(x, "sparseMatrix")) so the
comments in the above description of the coerce method
apply here too.
signature(a = "CHMfactor", b = "ddenseMatrix"), system= *:
The solve methods for a "CHMfactor" object take an
optional third argument system whose value can be one of the
character strings "A", "LDLt", "LD",
"DLt", "L", "Lt", "D", "P" or
"Pt". This argument describes the system to be solved. The
default, "A", is to solve \(Ax = b\) for \(x\) where
A is the sparse, positive-definite matrix that was factored
to produce a. Analogously, system = "L" returns the
solution \(x\), of \(Lx = b\).
Similarly, for all system codes but "P" and "Pt"
where, e.g., x <- solve(a, b, system="P") is equivalent to
x <- P %*% b.
See also solve-methods.
signature(x = "CHMfactor", logarithm =
"logical") returns the determinant (or the logarithm of the
determinant, if logarithm = TRUE, the default) of the
factor \(L\) from the \(LL'\) decomposition (even if the
decomposition represented by x is of the \(LDL'\)
form (!)). This is the square root of the determinant (half the
logarithm of the determinant when logarithm = TRUE) of the
positive-definite matrix that was decomposed.
%% since 0.999375-8 (2008-03-25):
signature(object = "CHMfactor"), parent. The
update method requires an additional argument
parent, which is either a
"dsCMatrix" object, say \(A\), (with the
same structure of nonzeros as the matrix that was decomposed to
produce object) or a general "dgCMatrix",
say \(M\), where \(A := M M'\) (== tcrossprod(parent))
is used for \(A\).
Further it provides an optional argument mult, a numeric
scalar. This method updates the numeric values in object
to the decomposition of \(A+mI\) where \(A\) is the matrix
above (either the parent or \(M M'\)) and \(m\) is
the scalar mult. Because only the numeric values are
updated this method should be faster than creating and decomposing
\(A+mI\). It is not uncommon to want, say, the determinant of
\(A+mI\) for many different values of \(m\). This method
would be the preferred approach in such cases.
Cholesky, also for examples;
class dgCMatrix.
## An example for the expand() method
n <- 1000; m <- 200; nnz <- 2000
set.seed(1)
M1 <- spMatrix(n, m,
i = sample(n, nnz, replace = TRUE),
j = sample(m, nnz, replace = TRUE),
x = round(rnorm(nnz),1))
XX <- crossprod(M1) ## = M1'M1 = M M' where M <- t(M1)
CX <- Cholesky(XX)
isLDL(CX)
str(CX) ## a "dCHMsimpl" object
r <- expand(CX)
L.P <- with(r, crossprod(L,P)) ## == L'P
PLLP <- crossprod(L.P) ## == (L'P)' L'P == P'LL'P = XX = M M'
b <- sample(m)
stopifnot(all.equal(PLLP, XX),
all(as.vector(solve(CX, b, system="P" )) == r$P %*% b),
all(as.vector(solve(CX, b, system="Pt")) == t(r$P) %*% b) )
u1 <- update(CX, XX, mult=pi)
u2 <- update(CX, t(M1), mult=pi) # with the original M, where XX = M M'
stopifnot(all.equal(u1,u2, tol=1e-14))
## [ See help(Cholesky) for more examples ]
## -------------
Run the code above in your browser using DataLab