Dense LU Factorizations

denseLU is the class of dense, row-pivoted LU factorizations of $m \times n$ real matrices $A$, having the general form $$P_{1} A = L U$$ or (equivalently) $$A = P_{1}' L U$$ where $P_{1}$ is an $m \times m$ permutation matrix, $L$ is an $m \times \min(m,n)$ unit lower trapezoidal matrix, and $U$ is a $\min(m,n) \times n$ upper trapezoidal matrix. If $m = n$, then the factors $L$ and $U$ are triangular.

Slots

Dim, Dimnames: inherited from virtual class MatrixFactorization.
x: a numeric vector of length prod(Dim) storing the triangular $L$ and $U$ factors together in a packed format. The details of the representation are specified by the manual for LAPACK routine dgetrf.
perm: an integer vector of length min(Dim) specifying the permutation $P_{1}$ as a product of transpositions. The corresponding permutation vector can be obtained as asPerm(perm).

Extends

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

Instantiation

Objects can be generated directly by calls of the form new("denseLU", ...), but they are more typically obtained as the value of lu(x) for x inheriting from denseMatrix (often dgeMatrix).

Methods

coerce: signature(from = "denseLU", to = "dgeMatrix"): returns a dgeMatrix with the dimensions of the factorized matrix $A$, equal to $L$ below the diagonal and equal to $U$ on and above the diagonal.
determinant: signature(from = "denseLU", logarithm = "logical"): computes the determinant of the factorized matrix $A$ or its logarithm.
expand: signature(x = "denseLU"): see expand-methods.
expand1: signature(x = "denseLU"): see expand1-methods.
expand2: signature(x = "denseLU"): see expand2-methods.
solve: signature(a = "denseLU", b = "missing"): see solve-methods.

References

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

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

Examples

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

n <- 3L
(A <- Matrix(round(rnorm(n * n), 2L), n, n))
#> 3 x 3 Matrix of class "dgeMatrix"
#>       [,1]  [,2] [,3]
#> [1,] -0.63  1.60 0.49
#> [2,]  0.18  0.33 0.74
#> [3,] -0.84 -0.82 0.58

## With dimnames, to see that they are propagated :
dimnames(A) <- dn <- list(paste0("r", seq_len(n)),
                          paste0("c", seq_len(n)))

(lu.A <- lu(A))
#> LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
#>   ..@ x       : num [1:9] -0.84 0.75 -0.214 -0.82 2.215 ...
#>   ..@ perm    : int [1:3] 3 3 3
#>   ..@ Dim     : int [1:2] 3 3
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:3] "r1" "r2" "r3"
#>   .. ..$ : chr [1:3] "c1" "c2" "c3"
str(e.lu.A <- expand2(lu.A), max.level = 2L)
#> List of 3
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#>  $ L  :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
#>  $ U  :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots

## Underlying LAPACK representation
(m.lu.A <- as(lu.A, "dgeMatrix")) # which is L and U interlaced
#> 3 x 3 Matrix of class "dgeMatrix"
#>            [,1]        [,2]      [,3]
#> [1,] -0.8400000 -0.82000000 0.5800000
#> [2,]  0.7500000  2.21500000 0.0550000
#> [3,] -0.2142857  0.06965495 0.8604547
stopifnot(identical(as(m.lu.A, "matrix"), `dim<-`(lu.A@x, lu.A@Dim)))

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

## A ~ P1' L U in floating point
stopifnot(exprs = {
    identical(names(e.lu.A), c("P1.", "L", "U"))
    identical(e.lu.A[["P1."]],
              new(  "pMatrix", Dim = c(n, n), Dimnames = c(dn[1L], list(NULL)),
                  margin = 1L, perm = invertPerm(asPerm(lu.A@perm))))
    identical(e.lu.A[["L"]],
              new("dtrMatrix", Dim = c(n, n), Dimnames = list(NULL, NULL),
                  uplo = "L", diag = "U", x = lu.A@x))
    identical(e.lu.A[["U"]],
              new("dtrMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
                  uplo = "U", diag = "N", x = lu.A@x))
    ae1(A, with(e.lu.A, P1. %*% L %*% U))
    ae2(A[asPerm(lu.A@perm), ], with(e.lu.A, L %*% U))
})

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