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.

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.

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

Examples

showClass("Cholesky")
#> Class "Cholesky" [package "Matrix"]
#> 
#> Slots:
#>                                                         
#> Name:       uplo         x      perm       Dim  Dimnames
#> Class: character   numeric   integer   integer      list
#> 
#> Extends: 
#> Class "CholeskyFactorization", directly
#> Class "MatrixFactorization", by class "CholeskyFactorization", distance 2
set.seed(1)

m <- 30L
n <- 6L
(A <- crossprod(Matrix(rnorm(m * n), m, n)))
#> 6 x 6 Matrix of class "dpoMatrix"
#>            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
#> [1,] 24.9699603  1.3658061  0.6673915  0.3077886  0.6145040 -7.6270511
#> [2,]  1.3658061 18.8724279 10.4969009 -2.0000204 -7.8836985  0.1955594
#> [3,]  0.6673915 10.4969009 27.1105459  0.3994417 -0.4068157  4.4884069
#> [4,]  0.3077886 -2.0000204  0.3994417 22.9601129 -3.7806013 -3.6029281
#> [5,]  0.6145040 -7.8836985 -0.4068157 -3.7806013 27.9716577 -9.8757902
#> [6,] -7.6270511  0.1955594  4.4884069 -3.6029281 -9.8757902 33.9927949

## 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
#> Cholesky factorization of Formal class 'Cholesky' [package "Matrix"] with 5 slots
#>   ..@ uplo    : chr "U"
#>   ..@ x       : num [1:36] 5.83 0 0 0 0 ...
#>   ..@ perm    : int [1:6] 6 3 5 1 4 2
#>   ..@ Dim     : int [1:2] 6 6
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "x1" "x2" "x3" "x4" ...
#>   .. ..$ : chr [1:6] "x1" "x2" "x3" "x4" ...
str(e.ch.A <- expand2(ch.A, LDL =  TRUE), max.level = 2L)
#> List of 5
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#>  $ L1 :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
#>  $ D  :Formal class 'ddiMatrix' [package "Matrix"] with 4 slots
#>  $ L1.:Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
#>  $ P1 :Formal class 'pMatrix' [package "Matrix"] with 4 slots
str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L)
#> List of 4
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#>  $ L  :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
#>  $ L. :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
#>  $ P1 :Formal class 'pMatrix' [package "Matrix"] with 4 slots

## Underlying LAPACK representation
(m.ch.A <- as(ch.A, "dtrMatrix")) # which is L', not L, because
#> 6 x 6 Matrix of class "dtrMatrix"
#>             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
#> [1,]  5.83033403  0.76983701 -1.69386353 -1.30816709 -0.61796255  0.03354172
#> [2,]           .  5.14955308  0.17422544  0.32516742  0.16995109  2.03339574
#> [3,]           .           .  5.00720776 -0.33112366 -0.96999278 -1.63387523
#> [4,]           .           .           .  4.80034190 -0.18270778  0.04322098
#> [5,]           .           .           .           .  4.64489868 -0.84002286
#> [6,]           .           .           .           .           .  3.37039314
A@uplo == "U"
#> [1] TRUE
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)
#> [1] TRUE
all.equal(solve(A, b), solve(ch.A, b), tolerance = 0)
#> [1] "Mean relative difference: 2.130541e-16"

## For identical results, we need the _unpivoted_ factorization
## computed by det(A) and solve(A, b)
(ch.A.nopivot <- Cholesky(A, perm = FALSE))
#> Cholesky factorization of Formal class 'Cholesky' [package "Matrix"] with 5 slots
#>   ..@ uplo    : chr "U"
#>   ..@ x       : num [1:36] 5 0 0 0 0 ...
#>   ..@ perm    : int(0) 
#>   ..@ Dim     : int [1:2] 6 6
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "x1" "x2" "x3" "x4" ...
#>   .. ..$ : chr [1:6] "x1" "x2" "x3" "x4" ...
stopifnot(identical(det(A), det(ch.A.nopivot)),
          identical(solve(A, b), solve(ch.A.nopivot, b)))