Cholesky-methods.RdComputes 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 non-negative and \(\sqrt{D}\) is well-defined.
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 class 'dsyMatrix'
Cholesky(A, perm = TRUE, tol = -1, ...)
# S4 method for class 'dspMatrix'
Cholesky(A, ...)
# S4 method for class 'dsCMatrix'
Cholesky(A, perm = TRUE, LDL = !super, super = FALSE,
Imult = 0, ...)
# S4 method for class 'ddiMatrix'
Cholesky(A, ...)
# S4 method for class 'generalMatrix'
Cholesky(A, uplo = "U", ...)
# S4 method for class 'triangularMatrix'
Cholesky(A, uplo = "U", ...)
# S4 method for class 'matrix'
Cholesky(A, uplo = "U", ...)a finite, symmetric matrix or
Matrix to be factorized. If A
is square but not symmetric, then it will be treated
as symmetric; see uplo.
Methods for dense A require positive definiteness
when perm = FALSE and positive semidefiniteness
when perm = TRUE.
Methods for sparse A require positive definiteness
when LDL = TRUE and nonzero leading principal minors
(after pivoting) when LDL = FALSE.
Methods for sparse, diagonal A are an exception,
requiring positive semidefiniteness unconditionally.
a logical indicating if the rows and columns
of \(A\) should be pivoted. Methods for sparse A
employ the approximate minimum degree (AMD) algorithm
in order to reduce fill-in, i.e., without regard for
numerical stability.
Pivoting for sparsity may introduce nonpositive leading
principal minors, causing the factorization to fail, in
which case it may be necessary to set perm = FALSE.
a finite numeric tolerance,
used only if perm = TRUE.
The factorization algorithm stops if the pivot is less than
or equal to tol. Negative tol is equivalent
to nrow(A) * .Machine$double.eps * max(diag(A)).
a logical indicating if the simplicial factorization
should be computed as
\(P_1' L_1 D L_1' P_1\),
such that the result stores the lower triangular entries
of \(L_1-I+D\).
The alternative is \(P_1' L L' P_1\),
such that the result stores the lower triangular entries
of \(L = L_1 \sqrt{D}\). This argument
is ignored if super = TRUE (or if super = NA
and the supernodal algorithm is chosen), as the supernodal
code does not yet support the LDL = TRUE variant.
a logical indicating if the factorization should
use the supernodal algorithm. The alternative is the simplicial
algorithm. Setting super = NA leaves the choice to
a CHOLMOD-internal heuristic.
a finite number. The matrix
that is factorized is A + Imult * diag(nrow(A)),
i.e., A plus Imult times the identity matrix.
This argument is useful for symmetric, indefinite A,
as Imult > max(rowSums(abs(A)) - diag(abs(A))) ensures
that A + Imult * diag(nrow(A)) is diagonally dominant.
(Symmetric, diagonally dominant matrices are positive definite.)
a string, either "U" or "L",
indicating which triangle of A should be used
to compute the factorization. The default is "U",
even for lower triangular A, to be consistent with
chol from base.
further arguments passed to or from methods.
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.
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.
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.
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). LAPACK-style 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, 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("Cholesky", inherited = FALSE)
#> Function: Cholesky (package Matrix)
#> A="ddiMatrix"
#> A="diagonalMatrix"
#> A="dsCMatrix"
#> A="dsRMatrix"
#> A="dsTMatrix"
#> A="dspMatrix"
#> A="dsyMatrix"
#> A="generalMatrix"
#> A="matrix"
#> A="symmetricMatrix"
#> A="triangularMatrix"
#>
set.seed(0)
## ---- Dense ----------------------------------------------------------
## .... Positive definite ..............................................
n <- 6L
(A1 <- crossprod(Matrix(rnorm(n * n), n, n)))
#> 6 x 6 Matrix of class "dpoMatrix"
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 7.63229784 3.5225440 -0.79842454 -0.04717246 -0.2763557 -2.99574643
#> [2,] 3.52254404 7.9529788 1.06842082 0.32822575 -2.7814087 -1.54583164
#> [3,] -0.79842454 1.0684208 2.51891543 -0.91363629 -0.4864619 -0.01953766
#> [4,] -0.04717246 0.3282258 -0.91363629 2.57854419 -1.2860880 1.44447913
#> [5,] -0.27635568 -2.7814087 -0.48646192 -1.28608796 3.5655733 -1.16142166
#> [6,] -2.99574643 -1.5458316 -0.01953766 1.44447913 -1.1614217 2.81492383
(ch.A1.nopivot <- Cholesky(A1, perm = FALSE))
#> Cholesky factorization of Formal class 'Cholesky' [package "Matrix"] with 5 slots
#> ..@ uplo : chr "U"
#> ..@ x : num [1:36] 2.76 0 0 0 0 ...
#> ..@ perm : int(0)
#> ..@ Dim : int [1:2] 6 6
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
(ch.A1 <- Cholesky(A1))
#> Cholesky factorization of Formal class 'Cholesky' [package "Matrix"] with 5 slots
#> ..@ uplo : chr "U"
#> ..@ x : num [1:36] 2.82 0 0 0 0 ...
#> ..@ perm : int [1:6] 2 1 4 5 3 6
#> ..@ Dim : int [1:2] 6 6
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
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)
#> 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.A1 <- expand2(ch.A1, 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
stopifnot(exprs = {
all.equal(as(A1, "matrix"), as(Reduce(`%*%`, e.ch.A1), "matrix"))
all.equal(as(A1, "matrix"), as(Reduce(`%*%`, E.ch.A1), "matrix"))
})
#> Error: as(A1, "matrix") and as(Reduce(`%*%`, e.ch.A1), "matrix") are not equal:
#> Mean relative difference: 1.275762
## .... Positive semidefinite but not positive definite ................
A2 <- A1
A2[1L, ] <- A2[, 1L] <- 0
A2
#> 6 x 6 Matrix of class "dgeMatrix"
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000
#> [2,] 0 7.9529788 1.06842082 0.3282258 -2.7814087 -1.54583164
#> [3,] 0 1.0684208 2.51891543 -0.9136363 -0.4864619 -0.01953766
#> [4,] 0 0.3282258 -0.91363629 2.5785442 -1.2860880 1.44447913
#> [5,] 0 -2.7814087 -0.48646192 -1.2860880 3.5655733 -1.16142166
#> [6,] 0 -1.5458316 -0.01953766 1.4444791 -1.1614217 2.81492383
try(Cholesky(A2, perm = FALSE)) # fails as not positive definite
#> Error in .local(A, ...) :
#> LAPACK routine 'dpotrf': leading principal minor of order 1 is not positive
ch.A2 <- Cholesky(A2) # returns, with a warning and ...
#> Warning: LAPACK routine 'dpstrf': matrix is rank deficient or not positive definite, the _computed_ rank is 5
A2.hat <- Reduce(`%*%`, expand2(ch.A2, LDL = FALSE))
norm(A2 - A2.hat, "2") / norm(A2, "2") # 7.670858e-17
#> [1] 0.6386362
## .... Not positive semidefinite ......................................
A3 <- A1
A3[1L, ] <- A3[, 1L] <- -1
A3
#> 6 x 6 Matrix of class "dgeMatrix"
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] -1 -1.0000000 -1.00000000 -1.0000000 -1.0000000 -1.00000000
#> [2,] -1 7.9529788 1.06842082 0.3282258 -2.7814087 -1.54583164
#> [3,] -1 1.0684208 2.51891543 -0.9136363 -0.4864619 -0.01953766
#> [4,] -1 0.3282258 -0.91363629 2.5785442 -1.2860880 1.44447913
#> [5,] -1 -2.7814087 -0.48646192 -1.2860880 3.5655733 -1.16142166
#> [6,] -1 -1.5458316 -0.01953766 1.4444791 -1.1614217 2.81492383
try(Cholesky(A3, perm = FALSE)) # fails as not positive definite
#> Error in .local(A, ...) :
#> LAPACK routine 'dpotrf': leading principal minor of order 1 is not positive
ch.A3 <- Cholesky(A3) # returns, with a warning and ...
#> Warning: LAPACK routine 'dpstrf': matrix is rank deficient or not positive definite, the _computed_ rank is 5
A3.hat <- Reduce(`%*%`, expand2(ch.A3, LDL = FALSE))
norm(A3 - A3.hat, "2") / norm(A3, "2") # 1.781568
#> [1] 0.7624035
## Indeed, 'A3' is not positive semidefinite, but 'A3.hat' _is_
ch.A3.hat <- Cholesky(A3.hat)
#> Warning: LAPACK routine 'dpstrf': matrix is rank deficient or not positive definite, the _computed_ rank is 3
A3.hat.hat <- Reduce(`%*%`, expand2(ch.A3.hat, LDL = FALSE))
norm(A3.hat - A3.hat.hat, "2") / norm(A3.hat, "2") # 1.777944e-16
#> [1] 0.91085
## ---- 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 )))
#> Warning: CHOLMOD warning 'matrix not positive definite' at file 'Supernodal/t_cholmod_super_numeric_worker.c', line 1114
#> Error in .local(A, ...): leading principal minor of order 81 is not positive
ch.A4
#> Error: object 'ch.A4' not found
s <- simplify2array
rapply2 <- function(object, f, ...) rapply(object, f, , , how = "list", ...)
s(rapply2(ch.A4, isLDL))
#> Error: object 'ch.A4' not found
s(m.ch.A4 <- rapply2(ch.A4, expand1, "L")) # giving L = L1 sqrt(D)
#> Error: object 'ch.A4' not found
## By design, the pivoted and simplicial factorizations
## are more sparse than the unpivoted and supernodal ones ...
s(rapply2(m.ch.A4, object.size))
#> Error: object 'm.ch.A4' not found
## Which is nicely visualized by lattice-based 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)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'print': error in evaluating the argument 'x' in selecting a method for function 'image': object 'm.ch.A4' not found
simpl1 <- ch.A4[[c("pivoted", "simpl1")]]
#> Error: object 'ch.A4' not found
stopifnot(exprs = {
length(simpl1@perm) == ncol(A4)
isPerm(simpl1@perm, 0L)
is.unsorted(simpl1@perm) # typically not the identity permutation
})
#> Error in eval(cl, envir = envir): object 'simpl1' not found
## 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)
#> Error: object 'simpl1' not found
stopifnot(min(diag(simpl1)) >= 0)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'diag': object 'simpl1' not found
str(e.ch.A4 <- expand2(simpl1, LDL = TRUE), max.level = 2L) # default
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'expand2': object 'simpl1' not found
str(E.ch.A4 <- expand2(simpl1, LDL = FALSE), max.level = 2L)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'expand2': object 'simpl1' not found
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"))
})
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'target' in selecting a method for function 'all.equal': object 'E.ch.A4' not found
## 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"]])
})
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'expand1': object 'simpl1' not found
## 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")
#> 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: 0x55c996d908c8>
#> <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"
stopifnot(all.equal(chol(A4, pivot = TRUE), E.ch.A4[["L."]]))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'current' in selecting a method for function 'all.equal': object 'E.ch.A4' not found
## 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)
#> [1] 1.398039e+02 3.058801e+00 2.680774e-01 -8.437545e-03 -1.396242e-01
#> [6] -3.356906e+00 -3.176258e+02
(t.ev <- table(factor(sign(ev), -1:1))) # the matrix "inertia"
#>
#> -1 0 1
#> 4 0 3
ch.A5 <- Cholesky(A5)
isLDL(ch.A5)
#> [1] TRUE
(d.A5 <- diag(ch.A5)) # diag(D) is partly negative
#> [1] -2 -4 -1 9 1 2 -1
## 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)
#> Error in expand1(ch.A5, "L") : D[i,i] is negative, i=1
try(expand2(ch.A5, LDL = FALSE)) # ditto
#> Error in .local(x, ...) : D[i,i] is negative, i=1
try(chol(A5, pivot = TRUE)) # ditto
#> Warning: CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449
#> Error in h(simpleError(msg, call)) :
#> error in evaluating the argument 'x' in selecting a method for function 't': leading principal minor of order 1 is not positive
## The default expansion is "square root free" and still works here
str(e.ch.A5 <- expand2(ch.A5, LDL = TRUE), max.level = 2L)
#> List of 5
#> $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#> $ L1 :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> $ D :Formal class 'ddiMatrix' [package "Matrix"] with 4 slots
#> $ L1.:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> $ P1 :Formal class 'pMatrix' [package "Matrix"] with 4 slots
stopifnot(all.equal(A5, as(Reduce(`%*%`, e.ch.A5), "symmetricMatrix")))
## Version of the SuiteSparse library, which includes CHOLMOD
Mv <- Matrix.Version()
Mv[["suitesparse"]]
#> [1] ‘7.6.0’