Description Usage Arguments Methods See Also Examples
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
1 2 | X <- solve(A,B)
|
In solve(a,b)
in the Matrix package, a
may also
be a MatrixFactorization
instead of directly a matrix.
1 2 3 4 5 6 7 8 9 | ## S4 method for signature 'CHMfactor,ddenseMatrix'
solve(a, b,
system = c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt"), ...)
## S4 method for signature '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
|
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 |
sparse |
only when |
tol |
only used when |
... |
potentially further arguments to the 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 argument system
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 where
A
is sparse, positive-definite matrix that was factored to produce
a
. 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 to x <- 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 via as(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
signature(a = "dgCMatrix", b = "ddenseMatrix")
with extra
argument list ( sparse = FALSE, tol = .Machine$double.eps )
:
Uses the sparse lu(a)
decomposition (which is cached
in a
's factor
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 if a
was sparse).
signature(a = "dgCMatrix", b = "dsparseMatrix")
, and
signature(a = "dgCMatrix", b = "missing")
with extra
argument list ( sparse=FALSE, tol = .Machine$double.eps )
:
Checks if a
is symmetric, and in that case, coerces it to
"symmetricMatrix"
, and then computes a
sparse solution via sparse Cholesky factorization,
independently of the sparse
argument. If a
is not
symmetric, the sparse lu
decomposition is used
and the result will be sparse or dense, depending on the
sparse
argument, exactly as for the above (b =
"ddenseMatrix"
) case.
signature(a = "dgeMatrix", b = ".....")
solve the system via internal LU, calling LAPACK routines
dgetri
or dgetrs
.
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
and dsytri
.
signature(a = "dtCMatrix", b = "CsparseMatrix")
,
signature(a = "dtCMatrix", b = "dgeMatrix")
, etc
sparse triangular solve, in traditional S/R also known as
backsolve
, or forwardsolve
.
solve(a,b)
is a sparseMatrix
if
b
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
, or forwardsolve
.
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")
.
solve
, lu
, and class documentations
CHMfactor
, sparseLU
, and
MatrixFactorization
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | ## 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)) )
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.