# Cholesky: Cholesky Decomposition of a Sparse Matrix In Matrix: Sparse and Dense Matrix Classes and Methods

## Description

Computes the Cholesky (aka “Choleski”) decomposition of a sparse, symmetric, positive-definite matrix. However, typically `chol()` should rather be used unless you are interested in the different kinds of sparse Cholesky decompositions.

## Usage

 `1` ```Cholesky(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...) ```

## Arguments

 `A` sparse symmetric matrix. No missing values or IEEE special values are allowed. `perm` logical scalar indicating if a fill-reducing permutation should be computed and applied to the rows and columns of `A`. Default is `TRUE`.
 `LDL` logical scalar indicating if the decomposition should be computed as LDL' where `L` is a unit lower triangular matrix. The alternative is LL' where `L` is lower triangular with arbitrary diagonal elements. Default is `TRUE`. Setting it to `NA` leaves the choice to a CHOLMOD-internal heuristic. `super` logical scalar indicating if a supernodal decomposition should be created. The alternative is a simplicial decomposition. Default is `FALSE`. Setting it to `NA` leaves the choice to a CHOLMOD-internal heuristic. `Imult` numeric scalar which defaults to zero. The matrix that is decomposed is A+m*I where m is the value of `Imult` and `I` is the identity matrix of order `ncol(A)`. `...` further arguments passed to or from other methods.

## Details

This is a generic function with special methods for different types of matrices. Use `showMethods("Cholesky")` to list all the methods for the `Cholesky` generic.

The method for class `dsCMatrix` of sparse matrices — the only one available currently — is based on functions from the CHOLMOD library.

Again: If you just want the Cholesky decomposition of a matrix in a straightforward way, you should probably rather use `chol(.)`.

Note that if `perm=TRUE` (default), the decomposition is

A = P' L~ D L~' P = P' L L' P,

where L can be extracted by `as(*, "Matrix")`, P by `as(*, "pMatrix")` and both by `expand(*)`, see the class `CHMfactor` documentation.

Note that consequently, you cannot easily get the “traditional” cholesky factor R, from this decomposition, as

R'R = A = P'LL'P = P' R~' R~ P = (R~ P)' (R~ P),

but R~ P is not triangular even though R~ is.

## Value

an object inheriting from either `"CHMsuper"`, or `"CHMsimpl"`, depending on the `super` argument; both classes extend `"CHMfactor"` which extends `"MatrixFactorization"`.

In other words, the result of `Cholesky()` is not a matrix, and if you want one, you should probably rather use `chol()`, see Details.

## References

Yanqing Chen, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam (2008) Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate. ACM Trans. Math. Softw. 35, 3, Article 22, 14 pages. doi: 10.1145/1391989.1391995

Timothy A. Davis (2006) Direct Methods for Sparse Linear Systems, SIAM Series “Fundamentals of Algorithms”.

Class definitions `CHMfactor` and `dsCMatrix` and function `expand`. Note the extra `solve(*, system = . )` options in `CHMfactor`.

Note that `chol()` returns matrices (inheriting from `"Matrix"`) whereas `Cholesky()` returns a `"CHMfactor"` object, and hence a typical user will rather use `chol(A)`.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58``` ```data(KNex) mtm <- with(KNex, crossprod(mm)) str(mtm@factors) # empty list() (C1 <- Cholesky(mtm)) # uses show() str(mtm@factors) # 'sPDCholesky' (simpl) (Cm <- Cholesky(mtm, super = TRUE)) c(C1 = isLDL(C1), Cm = isLDL(Cm)) str(mtm@factors) # 'sPDCholesky' *and* 'SPdCholesky' str(cm1 <- as(C1, "sparseMatrix")) str(cmat <- as(Cm, "sparseMatrix"))# hmm: super is *less* sparse here cm1[1:20, 1:20] b <- matrix(c(rep(0, 711), 1), nc = 1) ## solve(Cm, b) by default solves Ax = b, where A = Cm'Cm (= mtm)! ## hence, the identical() check *should* work, but fails on some GOTOblas: x <- solve(Cm, b) stopifnot(identical(x, solve(Cm, b, system = "A")), all.equal(x, solve(mtm, b))) Cn <- Cholesky(mtm, perm = FALSE)# no permutation -- much worse: sizes <- c(simple = object.size(C1), super = object.size(Cm), noPerm = object.size(Cn)) ## simple is 100, super= 137, noPerm= 812 : noquote(cbind(format(100 * sizes / sizes[1], digits=4))) ## Visualize the sparseness: dq <- function(ch) paste('"',ch,'"', sep="") ## dQuote() gives bad plots image(mtm, main=paste("crossprod(mm) : Sparse", dq(class(mtm)))) image(cm1, main= paste("as(Cholesky(crossprod(mm)),\"sparseMatrix\"):", dq(class(cm1)))) ## Smaller example, with same matrix as in help(chol) : (mm <- Matrix(toeplitz(c(10, 0, 1, 0, 3)), sparse = TRUE)) # 5 x 5 (opts <- expand.grid(perm = c(TRUE,FALSE), LDL = c(TRUE,FALSE), super = c(FALSE,TRUE))) rr <- lapply(seq_len(nrow(opts)), function(i) do.call(Cholesky, c(list(A = mm), opts[i,]))) nn <- do.call(expand.grid, c(attr(opts, "out.attr")\$dimnames, stringsAsFactors=FALSE,KEEP.OUT.ATTRS=FALSE)) names(rr) <- apply(nn, 1, function(r) paste(sub("(=.).*","\\1", r), collapse=",")) str(rr, max=1) str(re <- lapply(rr, expand), max=2) ## each has a 'P' and a 'L' matrix R0 <- chol(mm, pivot=FALSE) R1 <- chol(mm, pivot=TRUE ) stopifnot(all.equal(t(R1), re[[1]]\$L), all.equal(t(R0), re[[2]]\$L), identical(as(1:5, "pMatrix"), re[[2]]\$P), # no pivoting TRUE) # Version of the underlying SuiteSparse library by Tim Davis : .SuiteSparse_version() ```

### Example output

``` list()
'MatrixFactorization' of Formal class 'dCHMsimpl' [package "Matrix"] with 10 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 ...
..@ 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 ...
..@ type    : int [1:4] 2 0 0 1
..@ Dim     : int [1:2] 712 712
List of 1
\$ sPDCholesky:Formal class 'dCHMsimpl' [package "Matrix"] with 10 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 ...
.. ..@ 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 ...
.. ..@ type    : int [1:4] 2 0 0 1
.. ..@ Dim     : int [1:2] 712 712
'MatrixFactorization' of Formal class 'dCHMsuper' [package "Matrix"] with 9 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 ...
..@ 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 ...
..@ type    : int [1:6] 2 1 1 1 676 26
..@ Dim     : int [1:2] 712 712
C1    Cm
TRUE FALSE
List of 2
\$ sPDCholesky:Formal class 'dCHMsimpl' [package "Matrix"] with 10 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 ...
.. ..@ 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 ...
.. ..@ type    : int [1:4] 2 0 0 1
.. ..@ Dim     : int [1:2] 712 712
\$ SPdCholesky:Formal class 'dCHMsuper' [package "Matrix"] with 9 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 ...
.. ..@ 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 ...
.. ..@ type    : int [1:6] 2 1 1 1 676 26
.. ..@ Dim     : int [1:2] 712 712
Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
..@ i       : int [1:7451] 0 692 697 708 1 690 696 707 2 690 ...
..@ p       : int [1:713] 0 4 8 12 16 24 31 39 46 54 ...
..@ Dim     : int [1:2] 712 712
..@ Dimnames:List of 2
.. ..\$ : NULL
.. ..\$ : NULL
..@ x       : num [1:7451] 1 0.0919 -0.1241 -0.0984 1 ...
..@ uplo    : chr "L"
..@ diag    : chr "N"
Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
..@ i       : int [1:12501] 0 692 697 708 1 690 696 707 2 690 ...
..@ p       : int [1:713] 0 4 8 12 16 24 31 39 46 54 ...
..@ Dim     : int [1:2] 712 712
..@ Dimnames:List of 2
.. ..\$ : NULL
.. ..\$ : NULL
..@ x       : num [1:12501] 1 0.0919 -0.1241 -0.0984 1 ...
..@ uplo    : chr "L"
..@ diag    : chr "N"
20 x 20 sparse Matrix of class "dtCMatrix"

[1,] 1 . . . . . . . .            . . . . . . . . . . .
[2,] . 1 . . . . . . .            . . . . . . . . . . .
[3,] . . 1 . . . . . .            . . . . . . . . . . .
[4,] . . . 1 . . . . .            . . . . . . . . . . .
[5,] . . . . 1 . . . .            . . . . . . . . . . .
[6,] . . . . 0 1 . . .            . . . . . . . . . . .
[7,] . . . . . . 1 . .            . . . . . . . . . . .
[8,] . . . . . . 0 1 .            . . . . . . . . . . .
[9,] . . . . . . . . 1.000000e+00 . . . . . . . . . . .
[10,] . . . . . . . . 5.551115e-17 1 . . . . . . . . . .
[11,] . . . . . . . . .            . 1 . . . . . . . . .
[12,] . . . . . . . . .            . . 1 . . . . . . . .
[13,] . . . . . . . . .            . . . 1 . . . . . . .
[14,] . . . . . . . . .            . . . . 1 . . . . . .
[15,] . . . . . . . . .            . . . . . 1 . . . . .
[16,] . . . . . . . . .            . . . . . . 1 . . . .
[17,] . . . . . . . . .            . . . . . . . 1 . . .
[18,] . . . . . . . . .            . . . . . . . . 1 . .
[19,] . . . . . . . . .            . . . . . . . . . 1 .
[20,] . . . . . . . . .            . . . . . . . . . . 1
[,1]
simple 100.0
super  137.2
noPerm 811.9
5 x 5 sparse Matrix of class "dsCMatrix"

[1,] 10  .  1  .  3
[2,]  . 10  .  1  .
[3,]  1  . 10  .  1
[4,]  .  1  . 10  .
[5,]  3  .  1  . 10
perm   LDL super
1  TRUE  TRUE FALSE
2 FALSE  TRUE FALSE
3  TRUE FALSE FALSE
4 FALSE FALSE FALSE
5  TRUE  TRUE  TRUE
6 FALSE  TRUE  TRUE
7  TRUE FALSE  TRUE
8 FALSE FALSE  TRUE
List of 8
\$ perm=T,LDL=T,super=F:Formal class 'dCHMsimpl' [package "Matrix"] with 10 slots
\$ perm=F,LDL=T,super=F:Formal class 'dCHMsimpl' [package "Matrix"] with 10 slots
\$ perm=T,LDL=F,super=F:Formal class 'dCHMsimpl' [package "Matrix"] with 10 slots
\$ perm=F,LDL=F,super=F:Formal class 'dCHMsimpl' [package "Matrix"] with 10 slots
\$ perm=T,LDL=T,super=T:Formal class 'dCHMsuper' [package "Matrix"] with 9 slots
\$ perm=F,LDL=T,super=T:Formal class 'dCHMsuper' [package "Matrix"] with 9 slots
\$ perm=T,LDL=F,super=T:Formal class 'dCHMsuper' [package "Matrix"] with 9 slots
\$ perm=F,LDL=F,super=T:Formal class 'dCHMsuper' [package "Matrix"] with 9 slots
List of 8
\$ perm=T,LDL=T,super=F:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=F,LDL=T,super=F:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=T,LDL=F,super=F:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=F,LDL=F,super=F:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=T,LDL=T,super=T:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=F,LDL=T,super=T:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=T,LDL=F,super=T:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
\$ perm=F,LDL=F,super=T:List of 2
..\$ P:Formal class 'pMatrix' [package "Matrix"] with 4 slots
..\$ L:Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
[1] '4.2.1'
```

Matrix documentation built on June 11, 2021, 3 p.m.