Cholesky-class: Dense Cholesky Factorizations

Description

Classes Cholesky and pCholesky represent dense, pivoted Cholesky factorizations of $n \times n$ real, symmetric, positive semidefinite matrices $A$, having the general form $$P_{1} A P_{1}' = L_{1} D L_{1}' = L L'$$ or (equivalently) $$A = P_{1}' L_{1} D L_{1}' P_{1} = P_{1}' L L' P_{1}$$ where $P_{1}$ is a permutation matrix, $L_{1}$ is a unit lower triangular matrix, $D$ is a non-negative diagonal matrix, and $L = L_{1} \sqrt{D}$.

These classes store the entries of the Cholesky factor $L$ or its transpose $L'$ in a dense format as a vector of length $nn$ (Cholesky) or $n(n+1)/2$ (pCholesky), the latter giving the “packed” representation.

Arguments

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 $L'$ or $L$.
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 $L$ or its transpose $L'$ 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 $L$ or its transpose $L'$; see ‘Note’.
coerce: signature(from = "pCholesky", to = "dtpMatrix"): returns a dtpMatrix representing the Cholesky factor $L$ or its transpose $L'$; see ‘Note’.
determinant: signature(from = "p?Cholesky", logarithm = "logical"): computes the determinant of the factorized matrix $A$ or its logarithm.
diag: signature(x = "p?Cholesky"): returns a numeric vector of length $n$ containing the diagonal elements of $D$, which are the squared diagonal elements of $L$.
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.

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. tools:::Rd_expr_doi("10.56021/9781421407944")

Examples

Run this code

 
library(stats, pos = "package:base", verbose = FALSE)
library(utils, pos = "package:base", verbose = FALSE)

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

Run the code above in your browser using DataLab