| chol-methods {Matrix} | R Documentation |
Computes 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 signature 'dsyMatrix'
chol(x, pivot = FALSE, tol = -1, ...)
## S4 method for signature 'dspMatrix'
chol(x, ...)
## S4 method for signature 'dsCMatrix'
chol(x, pivot = FALSE, ...)
## S4 method for signature 'ddiMatrix'
chol(x, ...)
## S4 method for signature 'generalMatrix'
chol(x, uplo = "U", ...)
## S4 method for signature 'triangularMatrix'
chol(x, uplo = "U", ...)
x |
a finite, symmetric, positive
semidefinite matrix or |
pivot |
a logical indicating if the rows and columns
of |
tol |
a finite numeric tolerance,
used only if |
uplo |
a string, either |
... |
further arguments passed to or from methods. |
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}.
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.
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
The default method from base, chol,
called for traditional matrices x.
Generic function Cholesky, for more flexibility
notably when computing the Cholesky factorization and
not only the factor L'.
showMethods("chol", inherited = FALSE)
set.seed(0)
## ---- Dense ----------------------------------------------------------
## chol(x, pivot = value) wrapping Cholesky(x, perm = value)
selectMethod("chol", "dsyMatrix")
## Except in packed cases where pivoting is not yet available
selectMethod("chol", "dspMatrix")
## .... Positive definite ..............................................
(A1 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 5)))
(R1.nopivot <- chol(A1))
(R1 <- chol(A1, pivot = TRUE))
## 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)))
try(R2.nopivot <- chol(A2)) # fails as not positive definite
(R2 <- chol(A2, pivot = TRUE)) # returns, with a warning and ...
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)))
try(R3.nopivot <- chol(A3)) # fails as not positive definite
(R3 <- chol(A3, pivot = TRUE)) # returns, with a warning and ...
## _Not_ equal: see details and examples in help("Cholesky")
all.equal(t(A3[2:1, 2:1]), as(crossprod(R3), "dsyMatrix"))
## ---- Sparse ---------------------------------------------------------
## chol(x, pivot = value) wrapping
## Cholesky(x, perm = value, LDL = FALSE, super = FALSE)
selectMethod("chol", "dsCMatrix")
## Except in diagonal cases which are handled "directly"
selectMethod("chol", "ddiMatrix")
(A4 <- toeplitz(as(c(10, 0, 1, 0, 3), "sparseVector")))
(ch.A4.nopivot <- Cholesky(A4, perm = FALSE, LDL = FALSE, super = FALSE))
(ch.A4 <- Cholesky(A4, perm = TRUE, LDL = FALSE, super = FALSE))
(R4.nopivot <- chol(A4))
(R4 <- chol(A4, pivot = TRUE))
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))
})