Cholesky-class {Matrix}R Documentation

Dense Cholesky Factorizations

Description

Classes Cholesky and pCholesky represent dense, pivoted Cholesky factorizations of n×nn \times n real, symmetric, positive semidefinite matrices AA, having the general form

P1AP1=L1DL1=LLP_{1} A P_{1}' = L_{1} D L_{1}' = L L'

or (equivalently)

A=P1L1DL1P1=P1LLP1A = P_{1}' L_{1} D L_{1}' P_{1} = P_{1}' L L' P_{1}

where P1P_{1} is a permutation matrix, L1L_{1} is a unit lower triangular matrix, DD is a non-negative diagonal matrix, and L=L1DL = L_{1} \sqrt{D}.

These classes store the entries of the Cholesky factor LL or its transpose LL' in a dense format as a vector of length nnnn (Cholesky) or n(n+1)/2n(n+1)/2 (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 whether x stores LL' or LL.

x

a numeric vector of length n*n (Cholesky) or n*(n+1)/2 (pCholesky), where n=Dim[1], listing the entries of the Cholesky factor LL or its transpose LL' 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 a dtrMatrix representing the Cholesky factor LL or its transpose LL'; see ‘Note’.

coerce

signature(from = "pCholesky", to = "dtpMatrix"): returns a dtpMatrix representing the Cholesky factor LL or its transpose LL'; see ‘Note’.

determinant

signature(from = "p?Cholesky", logarithm = "logical"): computes the determinant of the factorized matrix AA or its logarithm.

diag

signature(x = "p?Cholesky"): returns a numeric vector of length nn containing the diagonal elements of DD, which are the squared diagonal elements of LL.

expand1

signature(x = "p?Cholesky"): see expand1-methods.

expand2

signature(x = "p?Cholesky"): see expand2-methods.

solve

signature(a = "p?Cholesky", b = .): see solve-methods.

Note

In Matrix < 1.6-0, class Cholesky extended dtrMatrix and class pCholesky extended dtpMatrix, reflecting the fact that the factor LL 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)))

[Package Matrix version 1.7-0 Index]