solve-methods: Methods in Package 'Matrix' for Function 'solve'
In Matrix: Sparse and Dense Matrix Classes and Methods

solve-methods

R Documentation

Methods in Package Matrix for Function `solve`

Description

Methods for generic function solve for solving linear systems of equations, i.e., for X in A X = B, where A is a square matrix and X and B are matrices with dimensions consistent with A.

Usage

solve(a, b, ...)

## S4 method for signature 'dgeMatrix,ANY'
solve(a, b, tol = .Machine$double.eps, ...)

## S4 method for signature 'dgCMatrix,missing'
solve(a, b, sparse = TRUE, ...)
## S4 method for signature 'dgCMatrix,matrix'
solve(a, b, sparse = FALSE, ...)
## S4 method for signature 'dgCMatrix,denseMatrix'
solve(a, b, sparse = FALSE, ...)
## S4 method for signature 'dgCMatrix,sparseMatrix'
solve(a, b, sparse = TRUE, ...)

## S4 method for signature 'denseLU,dgeMatrix'
solve(a, b, ...)
## S4 method for signature 'BunchKaufman,dgeMatrix'
solve(a, b, ...)
## S4 method for signature 'Cholesky,dgeMatrix'
solve(a, b, ...)
## S4 method for signature 'sparseLU,dgCMatrix'
solve(a, b, tol = .Machine$double.eps, ...)
## S4 method for signature 'sparseQR,dgCMatrix'
solve(a, b, ...)
## S4 method for signature 'CHMfactor,dgCMatrix'
solve(a, b, system = c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt"), ...)

Arguments

`a`	a finite square matrix or `Matrix` containing the coefficients of the linear system, or otherwise a `MatrixFactorization`, in which case methods behave (by default) as if the factorized matrix were specified.
`b`	a vector, `sparseVector`, matrix, or `Matrix` satisfying `NROW(b) == nrow(a)`, giving the right-hand side(s) of the linear system. Vectors `b` are treated as `length(b)`-by-1 matrices. If `b` is missing, then methods take `b` to be an identity matrix.
`tol`	a non-negative number. For `a` inheriting from `denseMatrix`, an error is signaled if the reciprocal one-norm condition number (see `rcond`) of `a` is less than `tol`, indicating that `a` is near-singular. For `a` of class `sparseLU`, an error is signaled if the ratio `min(d)/max(d)` is less than `tol`, where `d = abs(diag(a@U))`. (Interpret with care, as this ratio is a cheap heuristic and not in general equal to or even proportional to the reciprocal one-norm condition number.) Setting `tol = 0` disables the test.
`sparse`	a logical indicating if the result should be formally sparse, i.e., if the result should inherit from virtual class `sparseMatrix`. Only methods for sparse `a` and missing or matrix `b` have this argument. Methods for missing or sparse `b` use `sparse = TRUE` by default. Methods for dense `b` use `sparse = FALSE` by default.
`system`	a string specifying a linear system to be solved. Only methods for `a` inheriting from `CHMfactor` have this argument. See ‘Details’.
`...`	further arguments passed to or from methods.

Details

Methods for general and symmetric matrices a compute a triangular factorization (LU, Bunch-Kaufman, or Cholesky) and call the method for the corresponding factorization class. The factorization is sparse if a is. Methods for sparse, symmetric matrices a attempt a Cholesky factorization and perform an LU factorization only if that fails (typically because a is not positive definite).

Triangular, diagonal, and permutation matrices do not require factorization (they are already “factors”), hence methods for those are implemented directly. For triangular a, solutions are obtained by forward or backward substitution; for diagonal a, they are obtained by scaling the rows of b; and for permutations a, they are obtained by permuting the rows of b.

Methods for dense a are built on 14 LAPACK routines: class d..Matrix, where ..=(ge|tr|tp|sy|sp|po|pp), uses routines d..tri and d..trs for missing and non-missing b, respectively. A corollary is that these methods always give a dense result.

Methods for sparse a are built on CXSparse routines cs_lsolve, cs_usolve, and cs_spsolve and CHOLMOD routines cholmod_solve and cholmod_spsolve. By default, these methods give a vector result if b is a vector, a sparse matrix result if b is missing or a sparse matrix, and a dense matrix result if b is a dense matrix. One can override this behaviour by setting the sparse argument, where available, but that should be done with care. Note that a sparse result may be sparse only in the formal sense and not at all in the mathematical sense, depending on the nonzero patterns of a and b. Furthermore, whereas dense results are fully preallocated, sparse results must be “grown” in a loop over the columns of b.

Methods for a of class sparseQR are simple wrappers around qr.coef, giving the least squares solution in overdetermined cases.

Methods for a inheriting from CHMfactor can solve systems other than the default one A X = B. The correspondence between its system argument the system actually solved is outlined in the table below. See CHMfactor-class for a definition of notation.

`system`	`isLDL(a)=TRUE`	`isLDL(a)=FALSE`
`"A"`	`A X = B`	`A X = B`
`"LDLt"`	`L_{1} D L_{1}' X = B`	`L L' X = B`
`"LD"`	`L_{1} D X = B`	`L X = B`
`"DLt"`	`D L_{1}' X = B`	`L' X = B`
`"L"`	`L_{1} X = B`	`L X = B`
`"Lt"`	`L_{1}' X = B`	`L' X = B`
`"D"`	`D X = B`	`X = B`
`"P"`	`X = P_{1} B`	`X = P_{1} B`
`"Pt"`	`X = P_{1}' B`	`X = P_{1}' B`

Examples

## A close to symmetric example with "quite sparse" inverse:
n1 <- 7; n2 <- 3
dd <- data.frame(a = gl(n1,n2), b = gl(n2,1,n1*n2))# balanced 2-way
X <- sparse.model.matrix(~ -1+ a + b, dd)# no intercept --> even sparser
XXt <- tcrossprod(X)
diag(XXt) <- rep(c(0,0,1,0), length.out = nrow(XXt))

n <- nrow(ZZ <- kronecker(XXt, Diagonal(x=c(4,1))))
image(a <- 2*Diagonal(n) + ZZ %*% Diagonal(x=c(10, rep(1, n-1))))
isSymmetric(a) # FALSE
image(drop0(skewpart(a)))
image(ia0 <- solve(a, tol = 0)) # checker board, dense [but really, a is singular!]
try(solve(a, sparse=TRUE))##-> error [ TODO: assertError ]
ia. <- solve(a, sparse=TRUE, tol = 1e-19)##-> *no* error
if(R.version$arch == "x86_64")
  ## Fails on 32-bit [Fedora 19, R 3.0.2] from Matrix 1.1-0 on [FIXME ??] only
  stopifnot(all.equal(as.matrix(ia.), as.matrix(ia0)))
a <- a + Diagonal(n)
iad <- solve(a)
ias <- solve(a, sparse=FALSE)
stopifnot(all.equal(as(iad,"denseMatrix"), ias, tolerance=1e-14))
I. <- iad %*% a          ; image(I.)
I0 <- drop0(zapsmall(I.)); image(I0)
.I <- a %*% iad
.I0 <- drop0(zapsmall(.I))
stopifnot( all.equal(as(I0, "diagonalMatrix"), Diagonal(n)),
           all.equal(as(.I0,"diagonalMatrix"), Diagonal(n)) )

Matrix documentation built on Oct. 19, 2024, 1:08 a.m.