Methods for LU Factorization

Computes the pivoted LU factorization of an $m \times n$ real matrix $A$, which has the general form $$P_{1} A P_{2} = L U$$ or (equivalently) $$A = P_{1}' L U P_{2}'$$ where $P_{1}$ is an $m \times m$ permutation matrix, $P_{2}$ is an $n \times n$ 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.

Methods for denseMatrix are built on LAPACK routine dgetrf, which does not permute columns, so that $P_{2}$ is an identity matrix.

Methods for sparseMatrix are built on CXSparse routine cs_lu, which requires $m = n$, so that $L$ and $U$ are triangular matrices.

Usage

lu(x, ...)
# S4 method for class 'dgeMatrix'
lu(x, warnSing = TRUE, ...)
# S4 method for class 'dgCMatrix'
lu(x, errSing = TRUE, order = NA_integer_,
  tol = 1, ...)
# S4 method for class 'dsyMatrix'
lu(x, cache = TRUE, ...)
# S4 method for class 'dsCMatrix'
lu(x, cache = TRUE, ...)
# S4 method for class 'matrix'
lu(x, ...)

Arguments

x: a finite matrix or Matrix to be factorized, which must be square if sparse.
warnSing: a logical indicating if a warning should be signaled for singular x. Used only by methods for dense matrices.
errSing: a logical indicating if an error should be signaled for singular x. Used only by methods for sparse matrices.
order: an integer in 0:3 passed to CXSparse routine cs_sqr, indicating a strategy for choosing the column permutation $P_{2}$. 0 means no column permutation. 1, 2, and 3 indicate a fill-reducing ordering of $A + A'$, $\tilde{A}' \tilde{A}$, and $A' A$, where $\tilde{A}$ is $A$ with “dense” rows removed. NA (the default) is equivalent to 2 if tol == 1 and 1 otherwise. Do not set to 0 unless you know that the column order of $A$ is already sensible.
tol: a number. The original pivot element is used if its absolute value exceeds tol * a, where a is the maximum in absolute value of the other possible pivot elements. Set tol < 1 only if you know what you are doing.
cache: a logical indicating if the result should be cached in x@factors[["LU"]]. Note that caching is experimental and that only methods for classes extending generalMatrix or symmetricMatrix will have this argument.
...: further arguments passed to or from methods.

Value

An object representing the factorization, inheriting from virtual class LU. The specific class is denseLU unless x inherits from virtual class sparseMatrix, in which case it is sparseLU.

Details

What happens when x is determined to be near-singular differs by method. The method for class dgeMatrix completes the factorization, warning if warnSing = TRUE and in any case returning a valid denseLU object. Users of this method can detect singular x with a suitable warning handler; see tryCatch. In contrast, the method for class dgCMatrix abandons further computation, throwing an error if errSing = TRUE and otherwise returning NA. Users of this method can detect singular x with an error handler or by setting errSing = FALSE and testing for a formal result with is(., "sparseLU").

References

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

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

showMethods("lu", inherited = FALSE)
#> Function: lu (package Matrix)
#> x="denseMatrix"
#> x="dgCMatrix"
#> x="dgRMatrix"
#> x="dgTMatrix"
#> x="dgeMatrix"
#> x="diagonalMatrix"
#> x="dsCMatrix"
#> x="dsRMatrix"
#> x="dsTMatrix"
#> x="dspMatrix"
#> x="dsyMatrix"
#> x="dtCMatrix"
#> x="dtRMatrix"
#> x="dtTMatrix"
#> x="dtpMatrix"
#> x="dtrMatrix"
#> x="matrix"
#> x="sparseMatrix"
#> 
set.seed(0)

## ---- Dense ----------------------------------------------------------

(A1 <- Matrix(rnorm(9L), 3L, 3L))
#> 3 x 3 Matrix of class "dgeMatrix"
#>            [,1]       [,2]         [,3]
#> [1,]  1.2629543  1.2724293 -0.928567035
#> [2,] -0.3262334  0.4146414 -0.294720447
#> [3,]  1.3297993 -1.5399500 -0.005767173
(lu.A1 <- lu(A1))
#> LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
#>   ..@ x       : num [1:9] 1.33 0.95 -0.245 -1.54 2.735 ...
#>   ..@ perm    : int [1:3] 3 3 3
#>   ..@ Dim     : int [1:2] 3 3
#>   ..@ Dimnames:List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL

(A2 <- round(10 * A1[, -3L]))
#> 3 x 2 Matrix of class "dgeMatrix"
#>      [,1] [,2]
#> [1,]   13   13
#> [2,]   -3    4
#> [3,]   13  -15
(lu.A2 <- lu(A2))
#> LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
#>   ..@ x       : num [1:6] 13 1 -0.231 13 -28 ...
#>   ..@ perm    : int [1:2] 1 3
#>   ..@ Dim     : int [1:2] 3 2
#>   ..@ Dimnames:List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL

## A ~ P1' L U in floating point
str(e.lu.A2 <- expand2(lu.A2), max.level = 2L)
#> List of 3
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 slots
#>  $ L  :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
#>  $ U  :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
stopifnot(all.equal(A2, Reduce(`%*%`, e.lu.A2)))

## ---- Sparse ---------------------------------------------------------

A3 <- as(readMM(system.file("external/pores_1.mtx", package = "Matrix")),
         "CsparseMatrix")
(lu.A3 <- lu(A3))
#> 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
#>   .. ..$ : NULL
#>   .. ..$ : NULL

## A ~ P1' L U P2' in floating point
str(e.lu.A3 <- expand2(lu.A3), max.level = 2L)
#> List of 4
#>  $ P1.:Formal class 'pMatrix' [package "Matrix"] with 4 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 4 slots
stopifnot(all.equal(A3, Reduce(`%*%`, e.lu.A3)))