Learn R Programming

Matrix (version 1.5-1)

solve-methods: Methods in Package Matrix for Function solve()

Description

Methods for generic function solve, for solving linear systems of equations. These solve for \(X\) in $$A X = B$$ where \(A\) is a square matrix and \(X\) and \(B\) are matrices with compatible dimensions. The usual R syntax is


                        x <- solve(a, b, ...)
  

where b may also be a vector, in which case it is treated as a 1-column matrix. Methods support a inheriting from virtual classes Matrix and MatrixFactorization and b inheriting from virtual classes Matrix and sparseVector.

Usage

## solve(a, b, ...) # the two-argument version, almost always preferred to
## solve(a,    ...) # the *rarely needed* one-argument version

# S4 method for dgCMatrix,missing
solve(a, b, sparse = NA, ...)
# S4 method for dgCMatrix,matrix
solve(a, b, sparse = FALSE, ...)
# S4 method for dgCMatrix,denseMatrix
solve(a, b, sparse = FALSE, ...)
# S4 method for dgCMatrix,sparseMatrix
solve(a, b, sparse = NA, tol = .Machine$double.eps, ...)
# S4 method for CHMfactor,denseMatrix
solve(a, b, system = c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt"), ...)

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\).

sparse

only when a is a sparseMatrix: logical specifying if the result should also (formally) be sparse.

tol

only when a is a sparseMatrix and sparse is TRUE: an error is signaled if the ratio min(d)/max(d), where d = abs(diag(U)) and A = L U, is less than tol, indicating (near-)singular \(A\).

system

only when 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).

...

potentially further arguments to the methods.

Methods

signature(a = "ANY", b = "ANY")

is simply the base package's S3 generic solve.

%% This is copy-paste in CHMfactor-class.Rd {FIXME ?}

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, "generalMatrix"), 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 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 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

Run this code
## 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)) )

Run the code above in your browser using DataLab