Learn R Programming

Matrix (version 1.7-1)

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

See Also

Class CHMfactor for sparse Cholesky factorizations.

Classes dpoMatrix and dppMatrix.

Generic functions Cholesky, expand1 and expand2.

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