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
a
is 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
a
is a'>sparseMatrix
, i.e., typically a'>dgCMatrix
: logical specifying if the result should be a (formally) sparse matrix.- tol
only used when
a
is 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
solve
methods for a"'>CHMfactor"
object take an optional third argumentsystem
whose 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\) whereA
is 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
b
is a'>sparseMatrix
,system
is 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
'sfactor
slot). 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 ifa
was sparse).signature(a = "dgCMatrix", b = "dsparseMatrix")
, and
signature(a = "dgCMatrix", b = "missing")
with extra argument list
( sparse=FALSE, tol = .Machine$double.eps )
: Checks ifa
is symmetric, and in that case, coerces it to"'>symmetricMatrix"
, and then computes a sparse solution via sparse Cholesky factorization, independently of thesparse
argument. Ifa
is not symmetric, the sparselu
decomposition is used and the result will be sparse or dense, depending on thesparse
argument, exactly as for the above (b = "ddenseMatrix"
) case.signature(a = "dgeMatrix", b = ".....")
solve the system via internal LU, calling LAPACK routines
dgetri
ordgetrs
.signature(a = "diagonalMatrix", b = "matrix")
and other
b
s: 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
a
is 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
a
is 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
,dsytrs
anddsytri
.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'>sparseMatrix
ifb
is, and hence a'>denseMatrix
otherwise.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)) )
# }