solve-methods
Methods in Package Matrix for Function solve()
Methods for function solve to solve a linear system of
equations, or equivalently, solve for \(X\) in
$$A X = B$$
where \(A\) is a square matrix, and \(X\), \(B\) are matrices
or vectors (which are treated as 1-column matrices), and the R syntax
is
X <- solve(A,B)
In solve(a,b) in the Matrix package, a may also
be a '>MatrixFactorization instead of directly a matrix.
- Keywords
- methods
Usage
# S4 method for CHMfactor,ddenseMatrix
solve(a, b,
system = c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt"), …)# S4 method for dgCMatrix,matrix
solve(a, b, sparse = FALSE, tol = .Machine$double.eps, …)
solve(a, b, ...) ## *the* two-argument version, almost always preferred to
# solve(a) ## the *rarely* needed one-argument version
Arguments
- a
a square numeric matrix, \(A\), typically of one of the classes in Matrix. Logical matrices are coerced to corresponding numeric ones.
- b
numeric vector or matrix (dense or sparse) as RHS of the linear system \(Ax = b\).
- system
only if
ais a'>CHMfactor: character string indicating the kind of linear system to be solved, see below. Note that the default,"A", does not solve the triangular system (but"L"does).- sparse
only when
ais a'>sparseMatrix, i.e., typically a'>dgCMatrix: logical specifying if the result should be a (formally) sparse matrix.- tol
only used when
ais sparse, in theisSymmetric(a, tol=*)test, where that applies.- …
potentially further arguments to the methods.
Methods
signature(a = "ANY", b = "ANY")is simply the base package's S3 generic
solve.signature(a = "CHMfactor", b = "...."), system= *The
solvemethods for a"'>CHMfactor"object take an optional third argumentsystemwhose value can be one of the character strings"A","LDLt","LD","DLt","L","Lt","D","P"or"Pt". This argument describes the system to be solved. The default,"A", is to solve \(Ax = b\) for \(x\) whereAis sparse, positive-definite matrix that was factored to producea. Analogously,system = "L"returns the solution \(x\), of \(Lx = b\); similarly, for all system codes but"P"and"Pt"where, e.g.,x <- solve(a, b,system="P")is equivalent tox <- P %*% b.If
bis a'>sparseMatrix,systemis used as above the corresponding sparse CHOLMOD algorithm is called.signature(a = "ddenseMatrix", b = "....")(for all
b) work viaas(a, "dgeMatrix"), using the its methods, see below.signature(a = "denseLU", b = "missing")basically computes uses triangular forward- and back-solve.
signature(a = "dgCMatrix", b = "matrix"), and
%% -> ../R/dgCMatrix.R
signature(a = "dgCMatrix", b = "ddenseMatrix")with extra argument list
( sparse = FALSE, tol = .Machine$double.eps ): Uses the sparselu(a)decomposition (which is cached ina'sfactorslot). By default,sparse=FALSE, returns a'>denseMatrix, since \(U^{-1} L^{-1} B\) may not be sparse at all, even when \(L\) and \(U\) are.If
sparse=TRUE, returns a'>sparseMatrix(which may not be very sparse at all, even ifawas sparse).signature(a = "dgCMatrix", b = "dsparseMatrix"), and
signature(a = "dgCMatrix", b = "missing")with extra argument list
( sparse=FALSE, tol = .Machine$double.eps ): Checks ifais symmetric, and in that case, coerces it to"'>symmetricMatrix", and then computes a sparse solution via sparse Cholesky factorization, independently of thesparseargument. Ifais not symmetric, the sparseludecomposition is used and the result will be sparse or dense, depending on thesparseargument, exactly as for the above (b = "ddenseMatrix") case.signature(a = "dgeMatrix", b = ".....")solve the system via internal LU, calling LAPACK routines
dgetriordgetrs.signature(a = "diagonalMatrix", b = "matrix")and other
bs: Of course this is trivially implemented, as \(D^{-1}\) is diagonal with entries \(1 / D[i,i]\).signature(a = "dpoMatrix", b = "....Matrix"), and
signature(a = "dppMatrix", b = "....Matrix")The Cholesky decomposition of
ais calculated (if needed) while solving the system.signature(a = "dsCMatrix", b = "....")All these methods first try Cholmod's Cholesky factorization; if that works, i.e., typically if
ais positive semi-definite, it is made use of. Otherwise, the sparse LU decomposition is used as for the “general” matrices of class"dgCMatrix".signature(a = "dspMatrix", b = "...."), and
signature(a = "dsyMatrix", b = "....")all end up calling LAPACK routines
dsptri,dsptrs,dsytrsanddsytri.signature(a = "dtCMatrix", b = "CsparseMatrix"),
signature(a = "dtCMatrix", b = "dgeMatrix"), etc sparse triangular solve, in traditional S/R also known as
backsolve, orforwardsolve.solve(a,b)is a'>sparseMatrixifbis, and hence a'>denseMatrixotherwise.signature(a = "dtrMatrix", b = "ddenseMatrix"), and
signature(a = "dtpMatrix", b = "matrix"), and similar
b, including"missing", and"diagonalMatrix":all use LAPACK based versions of efficient triangular
backsolve, orforwardsolve.signature(a = "Matrix", b = "diagonalMatrix")works via
as(b, "CsparseMatrix").signature(a = "sparseQR", b = "ANY")simply uses
qr.coef(a, b).signature(a = "pMatrix", b = ".....")these methods typically use
crossprod(a,b), as the inverse of a permutation matrix is the same as its transpose.signature(a = "TsparseMatrix", b = "ANY")all work via
as(a, "CsparseMatrix").
%% This is copy-paste in CHMfactor-class.Rd {FIXME ?}
See Also
solve, lu, and class documentations
'>CHMfactor, '>sparseLU, and
'>MatrixFactorization.
Examples
# NOT RUN {
## 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)) # 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=TRUE)
stopifnot(all.equal(as(ias,"denseMatrix"), iad, 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)) )
# }