Learn R Programming

Matrix (version 1.7-4)

BunchKaufman-class: Dense Bunch-Kaufman Factorizations

Description

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.

Arguments

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. tools:::Rd_expr_doi("10.56021/9781421407944")

See Also

Class dsyMatrix and its packed counterpart.

Generic functions BunchKaufman, expand1, and expand2.

Examples

Run this code
 
library(stats, pos = "package:base", verbose = FALSE)
library(utils, pos = "package:base", verbose = FALSE)

showClass("BunchKaufman")
set.seed(1)

n <- 6L
(A <- forceSymmetric(Matrix(rnorm(n * n), n, n)))

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

(bk.A <- BunchKaufman(A))
str(e.bk.A <- expand2(bk.A, complete = FALSE), max.level = 2L)
str(E.bk.A <- expand2(bk.A, complete =  TRUE), max.level = 2L)

## Underlying LAPACK representation
(m.bk.A <- as(bk.A, "dtrMatrix"))
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)

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)))

Run the code above in your browser using DataLab