Choleskymethods  R Documentation 
Computes the pivoted Cholesky factorization of an
n \times n
real, symmetric matrix A
,
which has the general form
P_1 A P_1' = L_1 D L_1' \overset{D_{jj} \ge 0}{=} L L'
or (equivalently)
A = P_1' L_1 D L_1' P_1 \overset{D_{jj} \ge 0}{=} P_1' L L' P_1
where
P_1
is a permutation matrix,
L_1
is a unit lower triangular matrix,
D
is a diagonal matrix, and
L = L_1 \sqrt{D}
.
The second equalities hold only for positive semidefinite A
,
for which the diagonal entries of D
are nonnegative
and \sqrt{D}
is welldefined.
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
.
Cholesky(A, ...)
## S4 method for signature 'dsyMatrix'
Cholesky(A, perm = TRUE, tol = 1, ...)
## S4 method for signature 'dspMatrix'
Cholesky(A, ...)
## S4 method for signature 'dsCMatrix'
Cholesky(A, perm = TRUE, LDL = !super, super = FALSE,
Imult = 0, ...)
## S4 method for signature 'ddiMatrix'
Cholesky(A, ...)
## S4 method for signature 'generalMatrix'
Cholesky(A, uplo = "U", ...)
## S4 method for signature 'triangularMatrix'
Cholesky(A, uplo = "U", ...)
## S4 method for signature 'matrix'
Cholesky(A, uplo = "U", ...)
A 
a finite, symmetric matrix or

perm 
a logical indicating if the rows and columns
of 
tol 
a finite numeric tolerance,
used only if 
LDL 
a logical indicating if the simplicial factorization
should be computed as

super 
a logical indicating if the factorization should
use the supernodal algorithm. The alternative is the simplicial
algorithm. Setting 
Imult 
a finite number. The matrix
that is factorized is 
uplo 
a string, either 
... 
further arguments passed to or from methods. 
Note that the result of a call to Cholesky
inherits
from CholeskyFactorization
but not
Matrix
. Users who just want a matrix
should consider using chol
, whose methods are
simple wrappers around Cholesky
returning just the
upper triangular Cholesky factor L'
,
typically as a triangularMatrix
.
However, a more principled approach would be to construct
factors as needed from the CholeskyFactorization
object,
e.g., with expand1(x, "L")
, if x
is the
object.
The behaviour of Cholesky(A, perm = TRUE)
for dense A
is somewhat exceptional, in that it expects without checking
that A
is positive semidefinite. By construction, if A
is positive semidefinite and the exact algorithm encounters a zero
pivot, then the unfactorized trailing submatrix is the zero matrix,
and there is nothing left to do. Hence when the finite precision
algorithm encounters a pivot less than tol
, it signals a
warning instead of an error and zeros the trailing submatrix in
order to guarantee that P' L L' P
is positive
semidefinite even if A
is not. It follows that one way to
test for positive semidefiniteness of A
in the event of a
warning is to analyze the error
\frac{\lVert A  P' L L' P \rVert}{\lVert A \rVert}\,.
See the examples and LAPACK Working Note (“LAWN”) 161 for details.
An object representing the factorization, inheriting from
virtual class CholeskyFactorization
.
For a traditional matrix A
, the specific class is
Cholesky
.
For A
inheriting from
unpackedMatrix
,
packedMatrix
, and
sparseMatrix
,
the specific class is
Cholesky
,
pCholesky
, and
dCHMsimpl
or dCHMsuper
,
respectively.
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
.
Lucas, C. (2004). LAPACKstyle codes for level 2 and 3 pivoted Cholesky factorizations. LAPACK Working Note, Number 161. https://www.netlib.org/lapack/lawnspdf/lawn161.pdf
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, 114. \Sexpr[results=rd]{tools:::Rd_expr_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), 886905. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1145/1024074.1024081")}
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")}
Classes Cholesky
, pCholesky
,
dCHMsimpl
and dCHMsuper
and their methods.
Classes dpoMatrix
, dppMatrix
,
and dsCMatrix
.
Generic function chol
,
for obtaining the upper triangular Cholesky factor L'
as a
matrix or Matrix
.
Generic functions expand1
and expand2
,
for constructing matrix factors from the result.
Generic functions BunchKaufman
, Schur
,
lu
, and qr
,
for computing other factorizations.
showMethods("Cholesky", inherited = FALSE)
set.seed(0)
##  Dense 
## .... Positive definite ..............................................
n < 6L
(A1 < crossprod(Matrix(rnorm(n * n), n, n)))
(ch.A1.nopivot < Cholesky(A1, perm = FALSE))
(ch.A1 < Cholesky(A1))
stopifnot(exprs = {
length(ch.A1@perm) == ncol(A1)
isPerm(ch.A1@perm)
is.unsorted(ch.A1@perm) # typically not the identity permutation
length(ch.A1.nopivot@perm) == 0L
})
## A ~ P1' L D L' P1 ~ P1' L L' P1 in floating point
str(e.ch.A1 < expand2(ch.A1, LDL = TRUE), max.level = 2L)
str(E.ch.A1 < expand2(ch.A1, LDL = FALSE), max.level = 2L)
stopifnot(exprs = {
all.equal(as(A1, "matrix"), as(Reduce(`%*%`, e.ch.A1), "matrix"))
all.equal(as(A1, "matrix"), as(Reduce(`%*%`, E.ch.A1), "matrix"))
})
## .... Positive semidefinite but not positive definite ................
A2 < A1
A2[1L, ] < A2[, 1L] < 0
A2
try(Cholesky(A2, perm = FALSE)) # fails as not positive definite
ch.A2 < Cholesky(A2) # returns, with a warning and ...
A2.hat < Reduce(`%*%`, expand2(ch.A2, LDL = FALSE))
norm(A2  A2.hat, "2") / norm(A2, "2") # 7.670858e17
## .... Not positive semidefinite ......................................
A3 < A1
A3[1L, ] < A3[, 1L] < 1
A3
try(Cholesky(A3, perm = FALSE)) # fails as not positive definite
ch.A3 < Cholesky(A3) # returns, with a warning and ...
A3.hat < Reduce(`%*%`, expand2(ch.A3, LDL = FALSE))
norm(A3  A3.hat, "2") / norm(A3, "2") # 1.781568
## Indeed, 'A3' is not positive semidefinite, but 'A3.hat' _is_
ch.A3.hat < Cholesky(A3.hat)
A3.hat.hat < Reduce(`%*%`, expand2(ch.A3.hat, LDL = FALSE))
norm(A3.hat  A3.hat.hat, "2") / norm(A3.hat, "2") # 1.777944e16
##  Sparse 
## Really just three cases modulo permutation :
##
## type factorization minors of P1 A P1'
## 1 simplicial P1 A P1' = L1 D L1' nonzero
## 2 simplicial P1 A P1' = L L ' positive
## 3 supernodal P1 A P2' = L L ' positive
data(KNex, package = "Matrix")
A4 < crossprod(KNex[["mm"]])
ch.A4 <
list(pivoted =
list(simpl1 = Cholesky(A4, perm = TRUE, super = FALSE, LDL = TRUE),
simpl0 = Cholesky(A4, perm = TRUE, super = FALSE, LDL = FALSE),
super0 = Cholesky(A4, perm = TRUE, super = TRUE )),
unpivoted =
list(simpl1 = Cholesky(A4, perm = FALSE, super = FALSE, LDL = TRUE),
simpl0 = Cholesky(A4, perm = FALSE, super = FALSE, LDL = FALSE),
super0 = Cholesky(A4, perm = FALSE, super = TRUE )))
ch.A4
s < simplify2array
rapply2 < function(object, f, ...) rapply(object, f, , , how = "list", ...)
s(rapply2(ch.A4, isLDL))
s(m.ch.A4 < rapply2(ch.A4, expand1, "L")) # giving L = L1 sqrt(D)
## By design, the pivoted and simplicial factorizations
## are more sparse than the unpivoted and supernodal ones ...
s(rapply2(m.ch.A4, object.size))
## Which is nicely visualized by latticebased methods for 'image'
inm < c("pivoted", "unpivoted")
jnm < c("simpl1", "simpl0", "super0")
for(i in 1:2)
for(j in 1:3)
print(image(m.ch.A4[[c(i, j)]], main = paste(inm[i], jnm[j])),
split = c(j, i, 3L, 2L), more = i * j < 6L)
simpl1 < ch.A4[[c("pivoted", "simpl1")]]
stopifnot(exprs = {
length(simpl1@perm) == ncol(A4)
isPerm(simpl1@perm, 0L)
is.unsorted(simpl1@perm) # typically not the identity permutation
})
## One can expand with and without D regardless of isLDL(.),
## but "without" requires L = L1 sqrt(D), which is conditional
## on min(diag(D)) >= 0, hence "with" is the default
isLDL(simpl1)
stopifnot(min(diag(simpl1)) >= 0)
str(e.ch.A4 < expand2(simpl1, LDL = TRUE), max.level = 2L) # default
str(E.ch.A4 < expand2(simpl1, LDL = FALSE), max.level = 2L)
stopifnot(exprs = {
all.equal(E.ch.A4[["L" ]], e.ch.A4[["L1" ]] %*% sqrt(e.ch.A4[["D"]]))
all.equal(E.ch.A4[["L."]], sqrt(e.ch.A4[["D"]]) %*% e.ch.A4[["L1."]])
all.equal(A4, as(Reduce(`%*%`, e.ch.A4), "symmetricMatrix"))
all.equal(A4, as(Reduce(`%*%`, E.ch.A4), "symmetricMatrix"))
})
## The "same" permutation matrix with "alternate" representation
## [i, perm[i]] {margin=1} <> [invertPerm(perm)[j], j] {margin=2}
alt < function(P) {
P@margin < 1L + !(P@margin  1L) # 1 <> 2
P@perm < invertPerm(P@perm)
P
}
## Expansions are elegant but inefficient (transposes are redundant)
## hence programmers should consider methods for 'expand1' and 'diag'
stopifnot(exprs = {
identical(expand1(simpl1, "P1"), alt(e.ch.A4[["P1"]]))
identical(expand1(simpl1, "L"), E.ch.A4[["L"]])
identical(Diagonal(x = diag(simpl1)), e.ch.A4[["D"]])
})
## chol(A, pivot = value) is a simple wrapper around
## Cholesky(A, perm = value, LDL = FALSE, super = FALSE),
## returning L' = sqrt(D) L1' _but_ giving no information
## about the permutation P1
selectMethod("chol", "dsCMatrix")
stopifnot(all.equal(chol(A4, pivot = TRUE), E.ch.A4[["L."]]))
## Now a symmetric matrix with positive _and_ negative eigenvalues,
## hence _not_ positive semidefinite
A5 < new("dsCMatrix",
Dim = c(7L, 7L),
p = c(0:1, 3L, 6:7, 10:11, 15L),
i = c(0L, 0:1, 0:3, 2:5, 3:6),
x = c(1, 6, 38, 10, 60, 103, 4, 6, 32, 247, 2, 16, 128, 2, 67))
(ev < eigen(A5, only.values = TRUE)$values)
(t.ev < table(factor(sign(ev), 1:1))) # the matrix "inertia"
ch.A5 < Cholesky(A5)
isLDL(ch.A5)
(d.A5 < diag(ch.A5)) # diag(D) is partly negative
## Sylvester's law of inertia holds here, but not in general
## in finite precision arithmetic
stopifnot(identical(table(factor(sign(d.A5), 1:1)), t.ev))
try(expand1(ch.A5, "L")) # unable to compute L = L1 sqrt(D)
try(expand2(ch.A5, LDL = FALSE)) # ditto
try(chol(A5, pivot = TRUE)) # ditto
## The default expansion is "square root free" and still works here
str(e.ch.A5 < expand2(ch.A5, LDL = TRUE), max.level = 2L)
stopifnot(all.equal(A5, as(Reduce(`%*%`, e.ch.A5), "symmetricMatrix")))
## Version of the SuiteSparse library, which includes CHOLMOD
Mv < Matrix.Version()
Mv[["suitesparse"]]
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.