# 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 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

%% -> ../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

`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")`

.

%% 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)) )
# }
```

*Documentation reproduced from package Matrix, version 1.3-2, License: GPL (>= 2) | file LICENCE*