Cholesky-class: Dense Cholesky Factorizations

Cholesky-classR Documentation

Dense Cholesky Factorizations


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.


Dim, Dimnames

inherited from virtual class MatrixFactorization.


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.


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.


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.


Class CholeskyFactorization, directly. Class MatrixFactorization, by class CholeskyFactorization, distance 2.


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



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


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


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


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.


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


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


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


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.


The LAPACK source code, including documentation; see,, and

Lucas, C. (2004). LAPACK-style codes for level 2 and 3 pivoted Cholesky factorizations. LAPACK Working Note, Number 161.

Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.56021/9781421407944")}

See Also

Class CHMfactor for sparse Cholesky factorizations.

Classes dpoMatrix and dppMatrix.

Generic functions Cholesky, expand1 and expand2.



m <- 30L
n <- 6L
(A <- crossprod(Matrix(rnorm(m * n), m, n)))

## With dimnames, to see that they are propagated :
dimnames(A) <- dn <-"x", seq_len(n))), 2L)

(ch.A <- Cholesky(A)) # pivoted, by default
str( <- expand2(ch.A, LDL =  TRUE), max.level = 2L)
str( <- expand2(ch.A, LDL = FALSE), max.level = 2L)

## Underlying LAPACK representation
( <- as(ch.A, "dtrMatrix")) # which is L', not L, because
A@uplo == "U"
stopifnot(identical(as(, "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(, c("P1.", "L1", "D", "L1.", "P1"))
    identical(names(, c("P1.", "L" ,      "L." , "P1"))
              new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
                  margin = 2L, perm = invertPerm(ch.A@perm)))
    identical([["P1."]], t([["P1"]]))
    identical([["L1."]], t([["L1"]]))
    identical([["L." ]], t([["L" ]]))
    identical([["D"]], Diagonal(x = diag(ch.A)))
    all.equal([["L"]], with(, L1 %*% sqrt(D)))
    ae1(A, with(, P1. %*% L1 %*% D %*% L1. %*% P1))
    ae1(A, with(, P1. %*% L  %*%         L.  %*% P1))
    ae2(A[ch.A@perm, ch.A@perm], with(, L1 %*% D %*% L1.))
    ae2(A[ch.A@perm, ch.A@perm], with(, 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)))

Matrix documentation built on Nov. 14, 2023, 5:06 p.m.