R: Methods in Package Matrix for Function 'solve()'
solve-methods
R Documentation
Methods in Package Matrix for Function solve()
Description
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.
Usage
## 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
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 the
isSymmetric(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 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 awas 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 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 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").
See Also
solve, lu, and class documentations
CHMfactor, sparseLU, and
MatrixFactorization.
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)) # 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)) )
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Matrix)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Matrix/solve-methods.Rd_%03d_medium.png", width=480, height=480)
> ### Name: solve-methods
> ### Title: Methods in Package Matrix for Function 'solve()'
> ### Aliases: solve solve-methods solve,ANY,Matrix-method
> ### solve,CHMfactor,ANY-method solve,CHMfactor,ddenseMatrix-method
> ### solve,CHMfactor,diagonalMatrix-method
> ### solve,CHMfactor,dsparseMatrix-method solve,CHMfactor,matrix-method
> ### solve,CHMfactor,missing-method solve,CHMfactor,numeric-method
> ### solve,ddenseMatrix,ANY-method solve,ddenseMatrix,matrix-method
> ### solve,ddenseMatrix,Matrix-method solve,ddenseMatrix,missing-method
> ### solve,ddenseMatrix,numeric-method solve,denseLU,missing-method
> ### solve,dgCMatrix,ddenseMatrix-method
> ### solve,dgCMatrix,dsparseMatrix-method solve,dgCMatrix,matrix-method
> ### solve,dgCMatrix,missing-method solve,dgeMatrix,ddenseMatrix-method
> ### solve,dgeMatrix,matrix-method solve,dgeMatrix,missing-method
> ### solve,dgeMatrix,sparseMatrix-method
> ### solve,diagonalMatrix,matrix-method solve,diagonalMatrix,Matrix-method
> ### solve,diagonalMatrix,missing-method solve,dpoMatrix,dgeMatrix-method
> ### solve,dpoMatrix,matrix-method solve,dpoMatrix,missing-method
> ### solve,dppMatrix,dgeMatrix-method solve,dppMatrix,integer-method
> ### solve,dppMatrix,matrix-method solve,dppMatrix,missing-method
> ### solve,dsCMatrix,ddenseMatrix-method
> ### solve,dsCMatrix,denseMatrix-method
> ### solve,dsCMatrix,dsparseMatrix-method solve,dsCMatrix,matrix-method
> ### solve,dsCMatrix,missing-method solve,dsCMatrix,numeric-method
> ### solve,dspMatrix,ddenseMatrix-method solve,dspMatrix,matrix-method
> ### solve,dspMatrix,missing-method solve,dsyMatrix,ddenseMatrix-method
> ### solve,dsyMatrix,denseMatrix-method solve,dsyMatrix,matrix-method
> ### solve,dsyMatrix,missing-method solve,dtCMatrix,CsparseMatrix-method
> ### solve,dtCMatrix,dgeMatrix-method solve,dtCMatrix,matrix-method
> ### solve,dtCMatrix,missing-method solve,dtCMatrix,numeric-method
> ### solve,dtpMatrix,ddenseMatrix-method solve,dtpMatrix,matrix-method
> ### solve,dtpMatrix,missing-method solve,dtrMatrix,ddenseMatrix-method
> ### solve,dtrMatrix,dMatrix-method solve,dtrMatrix,matrix-method
> ### solve,dtrMatrix,Matrix-method solve,dtrMatrix,missing-method
> ### solve,Matrix,ANY-method solve,Matrix,diagonalMatrix-method
> ### solve,matrix,Matrix-method solve,Matrix,matrix-method
> ### solve,Matrix,missing-method solve,Matrix,numeric-method
> ### solve,Matrix,pMatrix-method solve,Matrix,sparseVector-method
> ### solve,MatrixFactorization,ANY-method
> ### solve,MatrixFactorization,missing-method
> ### solve,MatrixFactorization,numeric-method solve,pMatrix,matrix-method
> ### solve,pMatrix,Matrix-method solve,pMatrix,missing-method
> ### solve,sparseQR,ANY-method solve,TsparseMatrix,ANY-method
> ### solve,TsparseMatrix,missing-method
> ### Keywords: methods
>
> ### ** 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))))
Note: method with signature 'Matrix#diagonalMatrix' chosen for function 'kronecker',
target signature 'dsCMatrix#ddiMatrix'.
"sparseMatrix#ANY" would also be valid
Note: method with signature 'dsparseMatrix#dsparseMatrix' chosen for function 'kronecker',
target signature 'dsCMatrix#dtTMatrix'.
"sparseMatrix#TsparseMatrix" would also be valid
> image(a <- 2*Diagonal(n) + ZZ %*% Diagonal(x=c(10, rep(1, n-1))))
Note: method with signature 'TsparseMatrix#Matrix' chosen for function '%*%',
target signature 'dgTMatrix#ddiMatrix'.
"sparseMatrix#diagonalMatrix" would also be valid
> isSymmetric(a) # FALSE
[1] 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 ]
Error in .solve.dgC.lu(a, tol = tol) :
LU computationally singular: ratio of extreme entries in |diag(U)| = 3.036e-18
> 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)) )
>
>
>
>
>
>
> dev.off()
null device
1
>