Sparse LU Factorizations

sparseLU is the class of sparse, row- and column-pivoted LU factorizations of $n \times n$ real matrices $A$, having the general form $$P_{1} A P_{2} = L U$$ or (equivalently) $$A = P_{1}' L U P_{2}'$$ where $P_{1}$ and $P_{2}$ are permutation matrices, $L$ is a unit lower triangular matrix, and $U$ is an upper triangular matrix.

Slots

Dim, Dimnames: inherited from virtual class MatrixFactorization.
L: an object of class dtCMatrix, the unit lower triangular $L$ factor.
U: an object of class dtCMatrix, the upper triangular $U$ factor.
p, q: 0-based integer vectors of length Dim[1], specifying the permutations applied to the rows and columns of the factorized matrix. q of length 0 is valid and equivalent to the identity permutation, implying no column pivoting. Using R syntax, the matrix $P_{1} A P_{2}$ is precisely A[p+1, q+1] (A[p+1, ] when q has length 0).

Extends

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

Instantiation

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

Methods

determinant: signature(from = "sparseLU", logarithm = "logical"): computes the determinant of the factorized matrix $A$ or its logarithm.
expand: signature(x = "sparseLU"): see expand-methods.
expand1: signature(x = "sparseLU"): see expand1-methods.
expand2: signature(x = "sparseLU"): see expand2-methods.
solve: signature(a = "sparseLU", b = .): see solve-methods.

References

Davis, T. A. (2006). Direct methods for sparse linear systems. Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898718881

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

Examples

showClass("sparseLU")
#> Class "sparseLU" [package "Matrix"]
#> 
#> Slots:
#>                                                                   
#> Name:          L         U         p         q       Dim  Dimnames
#> Class: dtCMatrix dtCMatrix   integer   integer   integer      list
#> 
#> Extends: 
#> Class "LU", directly
#> Class "MatrixFactorization", by class "LU", distance 2
set.seed(2)

A <- as(readMM(system.file("external", "pores_1.mtx", package = "Matrix")),
        "CsparseMatrix")
(n <- A@Dim[1L])
#> [1] 30

## 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 'sparseLU' [package "Matrix"] with 6 slots
#>   ..@ L       :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>   .. .. ..@ i       : int [1:141] 0 1 18 22 25 26 1 18 22 25 ...
#>   .. .. ..@ p       : int [1:31] 0 6 11 13 16 19 22 26 31 35 ...
#>   .. .. ..@ Dim     : int [1:2] 30 30
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:141] 1 -0.993819 -0.004979 -0.000132 0.000132 ...
#>   .. .. ..@ uplo    : chr "L"
#>   .. .. ..@ diag    : chr "N"
#>   ..@ U       :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>   .. .. ..@ i       : int [1:206] 0 0 1 2 2 3 2 4 2 3 ...
#>   .. .. ..@ p       : int [1:31] 0 1 3 4 6 8 12 15 20 25 ...
#>   .. .. ..@ Dim     : int [1:2] 30 30
#>   .. .. ..@ Dimnames:List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : NULL
#>   .. .. ..@ x       : num [1:206] -7178502 -24613411 -18311731 -6399179 715 ...
#>   .. .. ..@ uplo    : chr "U"
#>   .. .. ..@ diag    : chr "N"
#>   ..@ p       : int [1:30] 1 11 29 27 19 28 9 18 17 7 ...
#>   ..@ q       : int [1:30] 0 1 29 27 19 28 9 18 17 8 ...
#>   ..@ Dim     : int [1:2] 30 30
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:30] "r1" "r2" "r3" "r4" ...
#>   .. ..$ : chr [1:30] "c1" "c2" "c3" "c4" ...
str(e.lu.A <- expand2(lu.A), max.level = 2L)
#> List of 4
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
#>  $ L  :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ U  :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#>  $ P2.:Formal class 'pMatrix' [package "Matrix"] with 5 slots

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 P2' in floating point
stopifnot(exprs = {
    identical(names(e.lu.A), c("P1.", "L", "U", "P2."))
    identical(e.lu.A[["P1."]],
              new("pMatrix", Dim = c(n, n), Dimnames = c(dn[1L], list(NULL)),
                  margin = 1L, perm = invertPerm(lu.A@p, 0L, 1L)))
    identical(e.lu.A[["P2."]],
              new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
                  margin = 2L, perm = invertPerm(lu.A@q, 0L, 1L)))
    identical(e.lu.A[["L"]], lu.A@L)
    identical(e.lu.A[["U"]], lu.A@U)
    ae1(A, with(e.lu.A, P1. %*% L %*% U %*% P2.))
    ae2(A[lu.A@p + 1L, lu.A@q + 1L], 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)))