chol.RdComputes the upper triangular Cholesky factor of an \(n \times n\) real, symmetric, positive semidefinite matrix \(A\), optionally after pivoting. That is the factor \(L'\) in $$P_{1} A P_{1}' = L L'$$ or (equivalently) $$A = P_{1}' L L' P_{1}$$ where \(P_{1}\) is a permutation matrix.
Methods for denseMatrix are built on
LAPACK routines dpstrf, dpotrf, and dpptrf,
The latter two do not permute rows or columns,
so that \(P_{1}\) is an identity matrix.
Methods for sparseMatrix are built on
CHOLMOD routines cholmod_analyze and cholmod_factorize_p.
chol(x, ...)
# S4 method for class 'dsyMatrix'
chol(x, pivot = FALSE, tol = -1, ...)
# S4 method for class 'dspMatrix'
chol(x, ...)
# S4 method for class 'dsCMatrix'
chol(x, pivot = FALSE, ...)
# S4 method for class 'ddiMatrix'
chol(x, ...)
# S4 method for class 'generalMatrix'
chol(x, uplo = "U", ...)
# S4 method for class 'triangularMatrix'
chol(x, uplo = "U", ...)a finite, symmetric, positive
semidefinite matrix or Matrix to
be factorized. If x is square but not symmetric,
then it will be treated as symmetric; see uplo.
Methods for dense x require positive definiteness
when pivot = FALSE.
Methods for sparse (but not diagonal) x require
positive definiteness unconditionally.
a logical indicating if the rows and columns
of \(x\) should be pivoted. Methods for sparse x
employ the approximate minimum degree (AMD) algorithm
in order to reduce fill-in, i.e., without regard for
numerical stability.
a finite numeric tolerance,
used only if pivot = TRUE.
The factorization algorithm stops if the pivot is less than
or equal to tol. Negative tol is equivalent
to nrow(x) * .Machine$double.eps * max(diag(x)).
a string, either "U" or "L",
indicating which triangle of x should be used
to compute the factorization. The default is "U",
even for lower triangular x, to be consistent with
chol from base.
further arguments passed to or from methods.
A matrix, triangularMatrix,
or diagonalMatrix representing
the upper triangular Cholesky factor \(L'\).
The result is a traditional matrix if x is a
traditional matrix, dense if x is dense, and
sparse if x is sparse.
For x inheriting from diagonalMatrix,
the diagonal result is computed directly and without pivoting,
i.e., bypassing CHOLMOD.
For all other x, chol(x, pivot = value) calls
Cholesky(x, perm = value, ...) under the hood.
If you must know the permutation \(P_{1}\) in addition
to the Cholesky factor \(L'\), then call Cholesky
directly, as the result of chol(x, pivot = TRUE) specifies
\(L'\) but not \(P_{1}\).
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.
The CHOLMOD source code; see
https://github.com/DrTimothyAldenDavis/SuiteSparse,
notably the header file CHOLMOD/Include/cholmod.h
defining cholmod_factor_struct.
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software, 35(3), Article 22, 1-14. doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004). Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software, 17(4), 886-905. doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
showMethods("chol", inherited = FALSE)
#> Function: chol (package base)
#> x="ANY"
#> x="ddiMatrix"
#> x="diagonalMatrix"
#> x="dsCMatrix"
#> x="dsRMatrix"
#> x="dsTMatrix"
#> x="dspMatrix"
#> x="dsyMatrix"
#> x="generalMatrix"
#> x="symmetricMatrix"
#> x="triangularMatrix"
#>
set.seed(0)
## ---- Dense ----------------------------------------------------------
## chol(x, pivot = value) wrapping Cholesky(x, perm = value)
selectMethod("chol", "dsyMatrix")
#> new("MethodDefinition", .Data = function (x, ...)
#> {
#> .local <- function (x, pivot = FALSE, tol = -1, ...)
#> {
#> ch <- as(Cholesky(x, perm = pivot, tol = tol), "dtrMatrix")
#> ch@Dimnames <- dimnames(x)
#> if (ch@uplo != "U")
#> t(ch)
#> else ch
#> }
#> .local(x, ...)
#> }, target = new("signature", .Data = "dsyMatrix", names = "x",
#> package = "Matrix"), defined = new("signature", .Data = "dsyMatrix",
#> names = "x", package = "Matrix"), generic = "chol")
#> <bytecode: 0x55b0434f4a10>
#> <environment: namespace:Matrix>
#> attr(,"target")
#> An object of class “signature”
#> x
#> "dsyMatrix"
#> attr(,"defined")
#> An object of class “signature”
#> x
#> "dsyMatrix"
#> attr(,"generic")
#> [1] "chol"
#> attr(,"generic")attr(,"package")
#> [1] "base"
#> attr(,"class")
#> [1] "MethodDefinition"
#> attr(,"class")attr(,"package")
#> [1] "methods"
## Except in packed cases where pivoting is not yet available
selectMethod("chol", "dspMatrix")
#> new("MethodDefinition", .Data = function (x, ...)
#> {
#> ch <- as(Cholesky(x), "dtpMatrix")
#> ch@Dimnames <- dimnames(x)
#> if (ch@uplo != "U")
#> t(ch)
#> else ch
#> }, target = new("signature", .Data = "dspMatrix", names = "x",
#> package = "Matrix"), defined = new("signature", .Data = "dspMatrix",
#> names = "x", package = "Matrix"), generic = "chol")
#> <bytecode: 0x55b0434f1738>
#> <environment: namespace:Matrix>
#> attr(,"target")
#> An object of class “signature”
#> x
#> "dspMatrix"
#> attr(,"defined")
#> An object of class “signature”
#> x
#> "dspMatrix"
#> attr(,"generic")
#> [1] "chol"
#> attr(,"generic")attr(,"package")
#> [1] "base"
#> attr(,"class")
#> [1] "MethodDefinition"
#> attr(,"class")attr(,"package")
#> [1] "methods"
## .... Positive definite ..............................................
(A1 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 5)))
#> 2 x 2 Matrix of class "dsyMatrix"
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 2 5
(R1.nopivot <- chol(A1))
#> 2 x 2 Matrix of class "dtrMatrix"
#> [,1] [,2]
#> [1,] 1 2
#> [2,] . 1
(R1 <- chol(A1, pivot = TRUE))
#> 2 x 2 Matrix of class "dtrMatrix"
#> [,1] [,2]
#> [1,] 2.2360680 0.8944272
#> [2,] . 0.4472136
## In 2-by-2 cases, we know that the permutation is 1:2 or 2:1,
## even if in general 'chol' does not say ...
stopifnot(exprs = {
all.equal( A1 , as(crossprod(R1.nopivot), "dsyMatrix"))
all.equal(t(A1[2:1, 2:1]), as(crossprod(R1 ), "dsyMatrix"))
identical(Cholesky(A1)@perm, 2:1) # because 5 > 1
})
## .... Positive semidefinite but not positive definite ................
(A2 <- new("dpoMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 4)))
#> 2 x 2 Matrix of class "dpoMatrix"
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 2 4
try(R2.nopivot <- chol(A2)) # fails as not positive definite
#> Error in .local(A, ...) :
#> LAPACK routine 'dpotrf': leading principal minor of order 2 is not positive
(R2 <- chol(A2, pivot = TRUE)) # returns, with a warning and ...
#> Warning: LAPACK routine 'dpstrf': matrix is rank deficient or not positive definite, the _computed_ rank is 1
#> 2 x 2 Matrix of class "dtrMatrix"
#> [,1] [,2]
#> [1,] 2 1
#> [2,] . 0
stopifnot(exprs = {
all.equal(t(A2[2:1, 2:1]), as(crossprod(R2), "dsyMatrix"))
identical(Cholesky(A2)@perm, 2:1) # because 4 > 1
})
## .... Not positive semidefinite ......................................
(A3 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 3)))
#> 2 x 2 Matrix of class "dsyMatrix"
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 2 3
try(R3.nopivot <- chol(A3)) # fails as not positive definite
#> Error in .local(A, ...) :
#> LAPACK routine 'dpotrf': leading principal minor of order 2 is not positive
(R3 <- chol(A3, pivot = TRUE)) # returns, with a warning and ...
#> Warning: LAPACK routine 'dpstrf': matrix is rank deficient or not positive definite, the _computed_ rank is 1
#> 2 x 2 Matrix of class "dtrMatrix"
#> [,1] [,2]
#> [1,] 1.732051 1.154701
#> [2,] . 0.000000
## _Not_ equal: see details and examples in help("Cholesky")
all.equal(t(A3[2:1, 2:1]), as(crossprod(R3), "dsyMatrix"))
#> [1] "Mean relative difference: 0.08333333"
## ---- Sparse ---------------------------------------------------------
## chol(x, pivot = value) wrapping
## Cholesky(x, perm = value, LDL = FALSE, super = FALSE)
selectMethod("chol", "dsCMatrix")
#> new("MethodDefinition", .Data = function (x, ...)
#> {
#> .local <- function (x, pivot = FALSE, ...)
#> {
#> ch <- t(as(Cholesky(x, perm = pivot, LDL = FALSE, super = FALSE),
#> "dtCMatrix"))
#> ch@Dimnames <- dimnames(x)
#> ch
#> }
#> .local(x, ...)
#> }, target = new("signature", .Data = "dsCMatrix", names = "x",
#> package = "Matrix"), defined = new("signature", .Data = "dsCMatrix",
#> names = "x", package = "Matrix"), generic = "chol")
#> <bytecode: 0x55b0434ecdd8>
#> <environment: namespace:Matrix>
#> attr(,"target")
#> An object of class “signature”
#> x
#> "dsCMatrix"
#> attr(,"defined")
#> An object of class “signature”
#> x
#> "dsCMatrix"
#> attr(,"generic")
#> [1] "chol"
#> attr(,"generic")attr(,"package")
#> [1] "base"
#> attr(,"class")
#> [1] "MethodDefinition"
#> attr(,"class")attr(,"package")
#> [1] "methods"
## Except in diagonal cases which are handled "directly"
selectMethod("chol", "ddiMatrix")
#> new("MethodDefinition", .Data = function (x, ...)
#> {
#> if (length(y <- x@x)) {
#> if (is.na(min.y <- min(y)) || min.y < 0)
#> stop(gettextf("%1$s(%2$s) is undefined: '%2$s' is not positive semidefinite",
#> "chol", "x"), domain = NA)
#> x@x <- sqrt(y)
#> }
#> x
#> }, target = new("signature", .Data = "ddiMatrix", names = "x",
#> package = "Matrix"), defined = new("signature", .Data = "ddiMatrix",
#> names = "x", package = "Matrix"), generic = "chol")
#> <bytecode: 0x55b0434ee508>
#> <environment: namespace:Matrix>
#> attr(,"target")
#> An object of class “signature”
#> x
#> "ddiMatrix"
#> attr(,"defined")
#> An object of class “signature”
#> x
#> "ddiMatrix"
#> attr(,"generic")
#> [1] "chol"
#> attr(,"generic")attr(,"package")
#> [1] "base"
#> attr(,"class")
#> [1] "MethodDefinition"
#> attr(,"class")attr(,"package")
#> [1] "methods"
(A4 <- toeplitz(as(c(10, 0, 1, 0, 3), "sparseVector")))
#> 5 x 5 sparse Matrix of class "dsCMatrix"
#>
#> [1,] 10 . 1 . 3
#> [2,] . 10 . 1 .
#> [3,] 1 . 10 . 1
#> [4,] . 1 . 10 .
#> [5,] 3 . 1 . 10
(ch.A4.nopivot <- Cholesky(A4, perm = FALSE, LDL = FALSE, super = FALSE))
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> ..@ x : num [1:9] 3.162 0.316 0.949 3.162 0.316 ...
#> ..@ p : int [1:6] 0 3 5 7 8 9
#> ..@ i : int [1:9] 0 2 4 1 3 2 4 3 4
#> ..@ nz : int [1:5] 3 2 2 1 1
#> ..@ nxt : int [1:7] 1 2 3 4 5 -1 0
#> ..@ prv : int [1:7] 6 0 1 2 3 4 -1
#> ..@ type : int [1:6] 0 1 0 1 0 0
#> ..@ colcount: int [1:5] 3 2 2 1 1
#> ..@ perm : int(0)
#> ..@ Dim : int [1:2] 5 5
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
(ch.A4 <- Cholesky(A4, perm = TRUE, LDL = FALSE, super = FALSE))
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> ..@ x : num [1:9] 3.162 0.316 3.146 3.162 0.316 ...
#> ..@ p : int [1:6] 0 2 3 6 8 9
#> ..@ i : int [1:9] 0 1 1 2 3 4 3 4 4
#> ..@ nz : int [1:5] 2 1 3 2 1
#> ..@ nxt : int [1:7] 1 2 3 4 5 -1 0
#> ..@ prv : int [1:7] 6 0 1 2 3 4 -1
#> ..@ type : int [1:6] 2 1 0 1 0 0
#> ..@ colcount: int [1:5] 2 1 3 2 1
#> ..@ perm : int [1:5] 1 3 0 2 4
#> ..@ Dim : int [1:2] 5 5
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
(R4.nopivot <- chol(A4))
#> 5 x 5 sparse Matrix of class "dtCMatrix"
#>
#> [1,] 3.162278 . 0.3162278 . 0.9486833
#> [2,] . 3.162278 . 0.3162278 .
#> [3,] . . 3.1464265 . 0.2224746
#> [4,] . . . 3.1464265 .
#> [5,] . . . . 3.0084057
(R4 <- chol(A4, pivot = TRUE))
#> 5 x 5 sparse Matrix of class "dtCMatrix"
#>
#> [1,] 3.162278 0.3162278 . . .
#> [2,] . 3.1464265 . . .
#> [3,] . . 3.162278 0.3162278 0.9486833
#> [4,] . . . 3.1464265 0.2224746
#> [5,] . . . . 3.0084057
det4 <- det(A4)
b4 <- rnorm(5L)
x4 <- solve(A4, b4)
stopifnot(exprs = {
identical(R4.nopivot, expand1(ch.A4.nopivot, "L."))
identical(R4, expand1(ch.A4, "L."))
all.equal(A4, crossprod(R4.nopivot))
all.equal(A4[ch.A4@perm + 1L, ch.A4@perm + 1L], crossprod(R4))
all.equal(diag(R4.nopivot), sqrt(diag(ch.A4.nopivot)))
all.equal(diag(R4), sqrt(diag(ch.A4)))
all.equal(sqrt(det4), det(R4.nopivot))
all.equal(sqrt(det4), det(R4))
all.equal(det4, det(ch.A4.nopivot, sqrt = FALSE))
all.equal(det4, det(ch.A4, sqrt = FALSE))
all.equal(x4, solve(R4.nopivot, solve(t(R4.nopivot), b4)))
all.equal(x4, solve(ch.A4.nopivot, b4))
all.equal(x4, solve(ch.A4, b4))
})