Cholesky-methods: Methods for Cholesky Factorization
In Matrix: Sparse and Dense Matrix Classes and Methods
Cholesky-methods R Documentation
Methods for Cholesky Factorization
Description
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 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.
Usage
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", ...)
Arguments
A
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.
perm
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.
tol
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)).
LDL
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.
super
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.
Imult
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.)
uplo
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.
Details
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.
Value
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.
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.
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.
\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), 886-905.
\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")}
See Also
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.
Examples
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.670858e-17
## .... 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.777944e-16
## ---- 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 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)
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"]]
Matrix documentation built on Oct. 19, 2024, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| Cholesky-methods | R Documentation |
Methods for Cholesky Factorization
Description
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 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.
Usage
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", ...)
Arguments
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. |
Details
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.
Value
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.
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.
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. \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), 886-905. \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")}
See Also
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.
Examples
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.670858e-17
## .... 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.777944e-16
## ---- 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 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)
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"]]
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.