Cholesky-class {Matrix} | R Documentation |
Dense Cholesky Factorizations
Description
Classes Cholesky
and pCholesky
represent
dense, pivoted Cholesky factorizations of
real, symmetric, positive semidefinite matrices
,
having the general form
or (equivalently)
where
is a permutation matrix,
is a unit lower triangular matrix,
is a non-negative diagonal matrix, and
.
These classes store the entries of the Cholesky factor
or its transpose
in a dense format as
a vector of length
(
Cholesky
) or
(
pCholesky
), the latter
giving the “packed” representation.
Slots
Dim
,Dimnames
inherited from virtual class
MatrixFactorization
.uplo
a string, either
"U"
or"L"
, indicating which triangle (upper or lower) of the factorized symmetric matrix was used to compute the factorization and in turn whetherx
storesor
.
x
a numeric vector of length
n*n
(Cholesky
) orn*(n+1)/2
(pCholesky
), wheren=Dim[1]
, listing the entries of the Cholesky factoror its transpose
in column-major order.
perm
a 1-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.
Extends
Class CholeskyFactorization
, directly.
Class MatrixFactorization
, by class
CholeskyFactorization
, distance 2.
Instantiation
Objects can be generated directly by calls of the form
new("Cholesky", ...)
or new("pCholesky", ...)
,
but they are more typically obtained as the value of
Cholesky(x)
for x
inheriting from
dsyMatrix
or dspMatrix
(often the subclasses of those reserved for positive
semidefinite matrices, namely dpoMatrix
and dppMatrix
).
Methods
coerce
signature(from = "Cholesky", to = "dtrMatrix")
: returns adtrMatrix
representing the Cholesky factoror its transpose
; see ‘Note’.
coerce
signature(from = "pCholesky", to = "dtpMatrix")
: returns adtpMatrix
representing the Cholesky factoror its transpose
; see ‘Note’.
determinant
signature(from = "p?Cholesky", logarithm = "logical")
: computes the determinant of the factorized matrixor its logarithm.
diag
signature(x = "p?Cholesky")
: returns a numeric vector of lengthcontaining the diagonal elements of
, which are the squared diagonal elements of
.
expand1
signature(x = "p?Cholesky")
: seeexpand1-methods
.expand2
signature(x = "p?Cholesky")
: seeexpand2-methods
.solve
signature(a = "p?Cholesky", b = .)
: seesolve-methods
.
Note
In Matrix < 1.6-0
, class Cholesky
extended
dtrMatrix
and class pCholesky
extended
dtpMatrix
, reflecting the fact that the factor
is indeed a triangular matrix.
Matrix
1.6-0
removed these extensions so that methods
would no longer be inherited from dtrMatrix
and dtpMatrix
.
The availability of such methods gave the wrong impression that
Cholesky
and pCholesky
represent a (singular)
matrix, when in fact they represent an ordered set of matrix factors.
The coercions as(., "dtrMatrix")
and as(., "dtpMatrix")
are provided for users who understand the caveats.
References
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dpstrf.f, https://netlib.org/lapack/double/dpotrf.f, and https://netlib.org/lapack/double/dpptrf.f.
Lucas, C. (2004). LAPACK-style codes for level 2 and 3 pivoted Cholesky factorizations. LAPACK Working Note, Number 161. https://www.netlib.org/lapack/lawnspdf/lawn161.pdf
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
See Also
Class CHMfactor
for sparse Cholesky factorizations.
Classes dpoMatrix
and dppMatrix
.
Generic functions Cholesky
,
expand1
and expand2
.
Examples
showClass("Cholesky")
set.seed(1)
m <- 30L
n <- 6L
(A <- crossprod(Matrix(rnorm(m * n), m, n)))
## 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)
## Underlying LAPACK representation
(m.ch.A <- as(ch.A, "dtrMatrix")) # which is L', not L, because
A@uplo == "U"
stopifnot(identical(as(m.ch.A, "matrix"), `dim<-`(ch.A@x, ch.A@Dim)))
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)))
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, ch.A@perm], with(e.ch.A, L1 %*% D %*% L1.))
ae2(A[ch.A@perm, ch.A@perm], with(E.ch.A, L %*% L. ))
})
## Factorization handled as factorized matrix
b <- rnorm(n)
all.equal(det(A), det(ch.A), tolerance = 0)
all.equal(solve(A, b), solve(ch.A, b), tolerance = 0)
## For identical results, we need the _unpivoted_ factorization
## computed by det(A) and solve(A, b)
(ch.A.nopivot <- Cholesky(A, perm = FALSE))
stopifnot(identical(det(A), det(ch.A.nopivot)),
identical(solve(A, b), solve(ch.A.nopivot, b)))