Cholesky-class function

Dense Cholesky Factorizations

Classes Cholesky and pCholesky represent dense, pivoted Cholesky factorizations of $n-by-n$

real, symmetric, positive semidefinite matrices $A$, having the general form [REMOVE_ME]$P_{1} A P_{1}' = L_{1} D L_{1}' = L L'P1 * A * P1' = L1 * D * L1' = L * L' [REMOVE_ME_2]$

or (equivalently) [REMOVE_ME]$A = P_{1}' L_{1} D L_{1}' P_{1} = P_{1}' L L' P_{1}A = P1' * L1 * D * L1' * P1 = P1' * L * L' * P1 [REMOVE_ME_2]$

where $P1$ is a permutation matrix, $L1$ is a unit lower triangular matrix, $D$ is a non-negative diagonal matrix, and $L = L1 * 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 $n*n$ (Cholesky) or $n*(n+1)/2$ (pCholesky), the latter giving the packed representation.

class

Description

Classes Cholesky and pCholesky represent dense, pivoted Cholesky factorizations of $n-by-n$

real, symmetric, positive semidefinite matrices $A$, having the general form

$$P_{1} A P_{1}' = L_{1} D L_{1}' = L L'P1 * A * P1' = L1 * D * L1' = L * L'$$

or (equivalently)

$$A = P_{1}' L_{1} D L_{1}' P_{1} = P_{1}' L L' P_{1}A = P1' * L1 * D * L1' * P1 = P1' * L * L' * P1$$

where $P1$ is a permutation matrix, $L1$ is a unit lower triangular matrix, $D$ is a non-negative diagonal matrix, and $L = L1 * 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 $n*n$ (Cholesky) or $n*(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 $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.

Note

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

See Also

Class CHMfactor for sparse Cholesky factorizations.

Classes dpoMatrix and dppMatrix.

Generic functions Cholesky, expand1 and expand2.

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

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