Methods for Cholesky Factorization
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.
Usage
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", ...)Arguments
- A
a finite, symmetric matrix or
Matrixto be factorized. IfAis square but not symmetric, then it will be treated as symmetric; seeuplo. Methods for denseArequire positive definiteness whenperm = FALSEand positive semidefiniteness whenperm = TRUE. Methods for sparseArequire positive definiteness whenLDL = TRUEand nonzero leading principal minors (after pivoting) whenLDL = FALSE. Methods for sparse, diagonalAare an exception, requiring positive semidefiniteness unconditionally.- perm
a logical indicating if the rows and columns of \(A\) should be pivoted. Methods for sparse
Aemploy 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 setperm = FALSE.- tol
a finite numeric tolerance, used only if
perm = TRUE. The factorization algorithm stops if the pivot is less than or equal totol. Negativetolis equivalent tonrow(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 ifsuper = NAand the supernodal algorithm is chosen), as the supernodal code does not yet support theLDL = TRUEvariant.- super
a logical indicating if the factorization should use the supernodal algorithm. The alternative is the simplicial algorithm. Setting
super = NAleaves the choice to a CHOLMOD-internal heuristic.- Imult
a finite number. The matrix that is factorized is
A + Imult * diag(nrow(A)), i.e.,AplusImulttimes the identity matrix. This argument is useful for symmetric, indefiniteA, asImult > max(rowSums(abs(A)) - diag(abs(A)))ensures thatA + Imult * diag(nrow(A))is diagonally dominant. (Symmetric, diagonally dominant matrices are positive definite.)- uplo
a string, either
"U"or"L", indicating which triangle ofAshould be used to compute the factorization. The default is"U", even for lower triangularA, to be consistent withcholfrom base.- ...
further arguments passed to or from methods.
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.
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.
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.
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. 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
Examples
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"))
})
## .... 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] 9.604062e-17
## .... 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.5916916
## 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 5
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] 4.841405e-17
## ---- 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
#> $pivoted
#> $pivoted$simpl1
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> ..@ x : num [1:7451] 1 0.0919 -0.1241 -0.0984 1 ...
#> ..@ p : int [1:713] 0 4 8 12 16 24 31 39 46 54 ...
#> ..@ i : int [1:7451] 0 692 697 708 1 690 696 707 2 690 ...
#> ..@ nz : int [1:712] 4 4 4 4 8 7 8 7 8 7 ...
#> ..@ nxt : int [1:714] 1 2 3 4 5 6 7 8 9 10 ...
#> ..@ prv : int [1:714] 713 0 1 2 3 4 5 6 7 8 ...
#> ..@ type : int [1:6] 2 0 0 1 0 0
#> ..@ colcount: int [1:712] 4 4 4 4 8 7 8 7 8 7 ...
#> ..@ perm : int [1:712] 256 243 242 241 213 693 212 692 125 633 ...
#> ..@ Dim : int [1:2] 712 712
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#>
#> $pivoted$simpl0
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> ..@ x : num [1:7451] 1 0.0919 -0.1241 -0.0984 1 ...
#> ..@ p : int [1:713] 0 4 8 12 16 24 31 39 46 54 ...
#> ..@ i : int [1:7451] 0 692 697 708 1 690 696 707 2 690 ...
#> ..@ nz : int [1:712] 4 4 4 4 8 7 8 7 8 7 ...
#> ..@ nxt : int [1:714] 1 2 3 4 5 6 7 8 9 10 ...
#> ..@ prv : int [1:714] 713 0 1 2 3 4 5 6 7 8 ...
#> ..@ type : int [1:6] 2 1 0 1 0 0
#> ..@ colcount: int [1:712] 4 4 4 4 8 7 8 7 8 7 ...
#> ..@ perm : int [1:712] 256 243 242 241 213 693 212 692 125 633 ...
#> ..@ Dim : int [1:2] 712 712
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#>
#> $pivoted$super0
#> Cholesky factorization of Formal class 'dCHMsuper' [package "Matrix"] with 10 slots
#> ..@ x : num [1:16616] 1 0.0919 -0.1241 -0.0984 1 ...
#> ..@ super : int [1:132] 0 1 2 3 4 6 8 10 11 12 ...
#> ..@ pi : int [1:132] 0 4 8 12 16 24 32 40 45 50 ...
#> ..@ px : int [1:132] 0 4 8 12 16 32 48 64 69 74 ...
#> ..@ s : int [1:1713] 0 692 697 708 1 690 696 707 2 690 ...
#> ..@ type : int [1:6] 2 1 1 1 676 26
#> ..@ colcount: int [1:712] 4 4 4 4 8 7 8 7 8 7 ...
#> ..@ perm : int [1:712] 256 243 242 241 213 693 212 692 125 633 ...
#> ..@ Dim : int [1:2] 712 712
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#>
#>
#> $unpivoted
#> $unpivoted$simpl1
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> ..@ x : num [1:71848] 1 0.2774 0.1961 0.2774 0.0958 ...
#> ..@ p : int [1:713] 0 16 22 36 55 63 76 82 90 101 ...
#> ..@ i : int [1:71848] 0 257 261 280 419 421 426 427 431 449 ...
#> ..@ nz : int [1:712] 16 6 14 19 8 13 6 8 11 9 ...
#> ..@ nxt : int [1:714] 1 2 3 4 5 6 7 8 9 10 ...
#> ..@ prv : int [1:714] 713 0 1 2 3 4 5 6 7 8 ...
#> ..@ type : int [1:6] 0 0 0 1 0 0
#> ..@ colcount: int [1:712] 16 6 14 19 8 13 6 8 11 9 ...
#> ..@ perm : int(0)
#> ..@ Dim : int [1:2] 712 712
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#>
#> $unpivoted$simpl0
#> Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> ..@ x : num [1:71848] 1 0.2774 0.1961 0.2774 0.0958 ...
#> ..@ p : int [1:713] 0 16 22 36 55 63 76 82 90 101 ...
#> ..@ i : int [1:71848] 0 257 261 280 419 421 426 427 431 449 ...
#> ..@ nz : int [1:712] 16 6 14 19 8 13 6 8 11 9 ...
#> ..@ nxt : int [1:714] 1 2 3 4 5 6 7 8 9 10 ...
#> ..@ prv : int [1:714] 713 0 1 2 3 4 5 6 7 8 ...
#> ..@ type : int [1:6] 0 1 0 1 0 0
#> ..@ colcount: int [1:712] 16 6 14 19 8 13 6 8 11 9 ...
#> ..@ perm : int(0)
#> ..@ Dim : int [1:2] 712 712
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#>
#> $unpivoted$super0
#> Cholesky factorization of Formal class 'dCHMsuper' [package "Matrix"] with 10 slots
#> ..@ x : num [1:139052] 1 0.2774 0.1961 0.2774 0.0958 ...
#> ..@ super : int [1:268] 0 1 2 3 4 5 6 7 8 9 ...
#> ..@ pi : int [1:268] 0 16 22 36 55 63 76 82 90 101 ...
#> ..@ px : int [1:268] 0 16 22 36 55 63 76 82 90 101 ...
#> ..@ s : int [1:3723] 0 257 261 280 419 421 426 427 431 449 ...
#> ..@ type : int [1:6] 0 1 1 1 44944 212
#> ..@ colcount: int [1:712] 16 6 14 19 8 13 6 8 11 9 ...
#> ..@ perm : int(0)
#> ..@ Dim : int [1:2] 712 712
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#>
#>
s <- simplify2array
rapply2 <- function(object, f, ...) rapply(object, f, , , how = "list", ...)
s(rapply2(ch.A4, isLDL))
#> pivoted unpivoted
#> simpl1 TRUE TRUE
#> simpl0 FALSE FALSE
#> super0 FALSE FALSE
s(m.ch.A4 <- rapply2(ch.A4, expand1, "L")) # giving L = L1 sqrt(D)
#> pivoted
#> simpl1 <S4 class ‘dtCMatrix’ [package “Matrix”] with 7 slots>
#> simpl0 <S4 class ‘dtCMatrix’ [package “Matrix”] with 7 slots>
#> super0 <S4 class ‘dgCMatrix’ [package “Matrix”] with 6 slots>
#> unpivoted
#> simpl1 <S4 class ‘dtCMatrix’ [package “Matrix”] with 7 slots>
#> simpl0 <S4 class ‘dtCMatrix’ [package “Matrix”] with 7 slots>
#> super0 <S4 class ‘dgCMatrix’ [package “Matrix”] with 6 slots>
## By design, the pivoted and simplicial factorizations
## are more sparse than the unpivoted and supernodal ones ...
s(rapply2(m.ch.A4, object.size))
#> pivoted unpivoted
#> simpl1 94056 bytes 866816 bytes
#> simpl0 94056 bytes 866816 bytes
#> super0 203744 bytes 1672976 bytes
## 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)
#> [1] TRUE
stopifnot(min(diag(simpl1)) >= 0)
str(e.ch.A4 <- expand2(simpl1, LDL = TRUE), max.level = 2L) # default
#> 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
str(E.ch.A4 <- expand2(simpl1, LDL = FALSE), max.level = 2L)
#> List of 4
#> $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#> $ L :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> $ L. :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> $ P1 :Formal class 'pMatrix' [package "Matrix"] with 4 slots
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")
#> 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: 0x564ca5020a70>
#> <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."]]))
## 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’