Dense Bunch-Kaufman Factorizations

Classes BunchKaufman and pBunchKaufman represent Bunch-Kaufman factorizations of $n \times n$ real, symmetric matrices $A$, having the general form $$A = U D_{U} U' = L D_{L} L'$$ where $D_{U}$ and $D_{L}$ are symmetric, block diagonal matrices composed of $b_{U}$ and $b_{L}$ $1 \times 1$ or $2 \times 2$ diagonal blocks; $U = \prod_{k = 1}^{b_{U}} P_{k} U_{k}$ is the product of $b_{U}$ row-permuted unit upper triangular matrices, each having nonzero entries above the diagonal in 1 or 2 columns; and $L = \prod_{k = 1}^{b_{L}} P_{k} L_{k}$ is the product of $b_{L}$ row-permuted unit lower triangular matrices, each having nonzero entries below the diagonal in 1 or 2 columns.

These classes store the nonzero entries of the $2 b_{U} + 1$ or $2 b_{L} + 1$ factors, which are individually sparse, in a dense format as a vector of length $nn$ (BunchKaufman) or $n(n+1)/2$ (pBunchKaufman), the latter giving the “packed” representation.

Slots

Dim, Dimnames: inherited from virtual class MatrixFactorization.
uplo: a string, either "U" or "L", indicating which triangle (upper or lower) of the factorized symmetric matrix was used to compute the factorization and in turn how the x slot is partitioned.
x: a numeric vector of length n*n (BunchKaufman) or n*(n+1)/2 (pBunchKaufman), where n=Dim[1]. The details of the representation are specified by the manual for LAPACK routines dsytrf and dsptrf.
perm: an integer vector of length n=Dim[1] specifying row and column interchanges as described in the manual for LAPACK routines dsytrf and dsptrf.

Extends

Class BunchKaufmanFactorization, directly. Class MatrixFactorization, by class BunchKaufmanFactorization, distance 2.

Instantiation

Objects can be generated directly by calls of the form new("BunchKaufman", ...) or new("pBunchKaufman", ...), but they are more typically obtained as the value of BunchKaufman(x) for x inheriting from dsyMatrix or dspMatrix.

Methods

coerce: signature(from = "BunchKaufman", to = "dtrMatrix"): returns a dtrMatrix, useful for inspecting the internal representation of the factorization; see ‘Note’.
coerce: signature(from = "pBunchKaufman", to = "dtpMatrix"): returns a dtpMatrix, useful for inspecting the internal representation of the factorization; see ‘Note’.
determinant: signature(from = "p?BunchKaufman", logarithm = "logical"): computes the determinant of the factorized matrix $A$ or its logarithm.
expand1: signature(x = "p?BunchKaufman"): see expand1-methods.
expand2: signature(x = "p?BunchKaufman"): see expand2-methods.
solve: signature(a = "p?BunchKaufman", b = .): see solve-methods.

Note

In Matrix < 1.6-0, class BunchKaufman extended dtrMatrix and class pBunchKaufman extended dtpMatrix, reflecting the fact that the internal representation of the factorization is fundamentally triangular: there are $n(n+1)/2$ “parameters”, and these can be arranged systematically to form an $n \times n$ triangular matrix. Matrix 1.6-0 removed these extensions so that methods would no longer be inherited from dtrMatrix and dtpMatrix. The availability of such methods gave the wrong impression that BunchKaufman and pBunchKaufman represent a (singular) matrix, when in fact they represent an ordered set of matrix factors.

The coercions as(., "dtrMatrix") and as(., "dtpMatrix") are provided for users who understand the caveats.

References

The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dsytrf.f and https://netlib.org/lapack/double/dsptrf.f.

Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944

Examples

showClass("BunchKaufman")
#> Class "BunchKaufman" [package "Matrix"]
#> 
#> Slots:
#>                                                         
#> Name:       uplo         x      perm       Dim  Dimnames
#> Class: character   numeric   integer   integer      list
#> 
#> Extends: 
#> Class "BunchKaufmanFactorization", directly
#> Class "MatrixFactorization", by class "BunchKaufmanFactorization", distance 2
set.seed(1)

n <- 6L
(A <- forceSymmetric(Matrix(rnorm(n * n), n, n)))
#> 6 x 6 Matrix of class "dsyMatrix"
#>            [,1]        [,2]       [,3]        [,4]        [,5]        [,6]
#> [1,] -0.6264538  0.48742905 -0.6212406  0.82122120  0.61982575  1.35867955
#> [2,]  0.4874291  0.73832471 -2.2146999  0.59390132 -0.05612874 -0.10278773
#> [3,] -0.6212406 -2.21469989  1.1249309  0.91897737 -0.15579551  0.38767161
#> [4,]  0.8212212  0.59390132  0.9189774  0.78213630 -1.47075238 -0.05380504
#> [5,]  0.6198257 -0.05612874 -0.1557955 -1.47075238 -0.47815006 -1.37705956
#> [6,]  1.3586796 -0.10278773  0.3876716 -0.05380504 -1.37705956 -0.41499456

## With dimnames, to see that they are propagated :
dimnames(A) <- rep.int(list(paste0("x", seq_len(n))), 2L)

(bk.A <- BunchKaufman(A))
#> Bunch-Kaufman factorization of Formal class 'BunchKaufman' [package "Matrix"] with 5 slots
#>   ..@ uplo    : chr "U"
#>   ..@ x       : num [1:36] -0.726 0 0 0 0 ...
#>   ..@ perm    : int [1:6] 1 -2 -2 4 -5 -5
#>   ..@ Dim     : int [1:2] 6 6
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:6] "x1" "x2" "x3" "x4" ...
#>   .. ..$ : chr [1:6] "x1" "x2" "x3" "x4" ...
str(e.bk.A <- expand2(bk.A, complete = FALSE), max.level = 2L)
#> List of 3
#>  $ U :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>  $ DU:Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
#>  $ U.:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
str(E.bk.A <- expand2(bk.A, complete =  TRUE), max.level = 2L)
#> List of 17
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ :Formal class 'pMatrix' [package "Matrix"] with 5 slots

## Underlying LAPACK representation
(m.bk.A <- as(bk.A, "dtrMatrix"))
#> 6 x 6 Matrix of class "dtrMatrix"
#>      [,1]        [,2]        [,3]        [,4]        [,5]        [,6]       
#> [1,] -0.72555450 -0.05748438 -0.49161210 -1.53065172 -0.95046421 -0.12008285
#> [2,]           . -0.53207259 -2.96649303  1.83012123  0.06964736  0.01657652
#> [3,]           .           .  0.53053617  1.08460713 -0.35250257  0.23553420
#> [4,]           .           .           .  0.38097326 -0.31584491  1.17770770
#> [5,]           .           .           .           . -0.47815006 -1.37705956
#> [6,]           .           .           .           .           . -0.41499456
stopifnot(identical(as(m.bk.A, "matrix"), `dim<-`(bk.A@x, bk.A@Dim)))

## Number of factors is 2*b+1, b <= n, which can be nontrivial ...
(b <- (length(E.bk.A) - 1L) %/% 2L)
#> [1] 8

ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...)
ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...)

## A ~ U DU U', U := prod(Pk Uk) in floating point
stopifnot(exprs = {
    identical(names(e.bk.A), c("U", "DU", "U."))
    identical(e.bk.A[["U" ]], Reduce(`%*%`, E.bk.A[seq_len(b)]))
    identical(e.bk.A[["U."]], t(e.bk.A[["U"]]))
    ae1(A, with(e.bk.A, U %*% DU %*% U.))
})

## Factorization handled as factorized matrix
b <- rnorm(n)
stopifnot(identical(det(A), det(bk.A)),
          identical(solve(A, b), solve(bk.A, b)))