SparseM.solve {SparseM}

Linear Equation Solving via Cholesky Decomposition for Sparse Matrices

Description

chol()

performs a Cholesky decomposition of a symmetric positive definite sparse matrix x of class matrix.csr.

backsolve()

performs a triangular back-fitting to compute the solutions of a system of linear equations in one step.

backsolve() and forwardsolve()

can also split the functionality of backsolve into two steps.

solve()

combines chol() and backsolve() to compute the inverse of a matrix if the right-hand-side is missing.

Usage

chol(x, ...)
## S4 method for signature 'matrix.csr'
chol(x, pivot = FALSE,
    nsubmax, nnzlmax, tmpmax,
    eps = .Machine$double.eps, tiny = 1e-30, Large = 1e128, warnOnly = FALSE,
    cacheKb = 1024L, level = 8L, ...)

## S4 method for signature 'matrix.csr.chol'
backsolve(r, x, k, upper.tri, transpose,
          twice = TRUE, drop = TRUE, ...)
## S4 method for signature 'matrix.csr.chol'
forwardsolve(l, x, k, upper.tri, transpose)
## S4 method for signature 'matrix.csr'
solve(a, b, ...)

Arguments

a

symmetric positive definite matrix of class "matrix.csr".

r, l

object of class "matrix.csr.chol" as returned by the chol() method.

x

For chol():

One of the sparse matrix classes, "matrix.csr" or "matrix.csc";

For {back,forward,}solve():

vector or regular matrix of right-hand-side(s) of a system of linear equations.

b

vector or matrix right-hand-side(s) to solve for.

k

inherited from the generic; not used here.

pivot

inherited from the generic; not used here.

nsubmax, nnzlmax, tmpmax

positive integer numbers with smart defaults; do not set unless you know what you are doing!

eps

positive tolerance for checking symmetry; change with caution.

tiny

positive tolerance for checking diagonal entries to be “essentially zero” and hence to be replaced by Large, during Cholesky decomposition. Chaning this value may help in close to singular cases, see ‘Examples’.

Large

large positive number, “essentially infinite”, to replace tiny diagonal entries during Cholesky.

warnOnly

logical; when set to true, a result is returned with a warning instead of an error (via stop()); notably in close to singular cases.

cacheKb

a positive integer, specifying an approximate size of the machine's cache memory in kilo (1024) bytes (‘Kb’); used to be hard wired to 64.

level

level of loop unrolling while performing numerical factorization; an integer in c(1, 2, 4, 8); used to be hard wired to 8.

upper.tri, transpose

inherited from the generic; not used here.

twice

logical flag: If true, backsolve() solves twice, see below.

drop

logical flag: If true, backsolve() returns drop(.), i.e., a vector instead of a column-1 matrix.

...

further arguments passed to or from other methods.

Details

chol performs a Cholesky decomposition of a symmetric positive definite sparse matrix a of class matrix.csr using the block sparse Cholesky algorithm of Ng and Peyton (1993). The structure of the resulting matrix.csr.chol object is relatively complicated. If necessary it can be coerced back to a matrix.csr object as usual with as.matrix.csr. backsolve does triangular back-fitting to compute the solutions of a system of linear equations. For systems of linear equations that only vary on the right-hand-side, the result from chol can be reused. Contrary to the behavior of backsolve in base R, the default behavior of backsolve(C,b) when C is a matrix.csr.chol object is to produce a solution to the system Ax = b where C <- chol(A), see the example section. When the flag twice is FALSE then backsolve solves the system Cx = b, up to a permutation – see the comments below. The command solve combines chol and backsolve, and will compute the inverse of a matrix if the right-hand-side is missing. The determinant of the Cholesky factor is returned providing a means to efficiently compute the determinant of sparse positive definite symmetric matrices.

There are several integer storage parameters that are set by default in the call to the Cholesky factorization, these can be overridden in any of the above functions and will be passed by the usual "dots" mechanism. The necessity to do this is usually apparent from error messages like: Error in local(X...) increase tmpmax. For example, one can use, solve(A,b, tmpmax = 100*nrow(A)). The current default for tmpmax is 50*nrow(A). Some experimentation may be needed to select appropriate values, since they are highly problem dependent. See the code of chol() for further details on the current defaults.

Note

There is no explicit checking for positive definiteness of the matrix so users are advised to ensure that this condition is satisfied. Messages such as "insufficient space" may indicate that one is trying to factor a singular matrix. Because the sparse Cholesky algorithm re-orders the positive definite sparse matrix A, the value of x <- backsolve(C, b) does not equal the solution to the triangular system Cx = b, but is instead the solution to the system CPx = Pb for some permutation matrix P (and analogously for x <- forwardsolve(C, b)). However, a little algebra easily shows that backsolve(C, forwardsolve(C, b), twice = FALSE) is the solution to the equation Ax=b. Finally, if C <- chol(A) for some sparse covariance matrix A, and z is a conformable standard normal vector, then the product y <- as.matrix.csr(C) %*% z is normal with covariance matrix A irrespective of the permutation of the Cholesky factor.

References

Koenker, R and Ng, P. (2002) SparseM: A Sparse Matrix Package for R. http://www.econ.uiuc.edu/~roger/research/home.html

Ng, E. G. and B. W. Peyton (1993) Block sparse Cholesky algorithms on advanced uniprocessor computers. SIAM J. Scientific Computing 14, 1034–1056.

See Also

slm() for a sparse version of stats package's lm().

Examples

data(lsq)
class(lsq) # -> [1] "matrix.csc.hb"
model.matrix(lsq)->design.o
class(design.o) # -> "matrix.csr"
dim(design.o) # -> [1] 1850  712
y <- model.response(lsq) # extract the rhs
length(y) # [1] 1850

X <- as.matrix(design.o)
c(object.size(X) / object.size(design.o)) ## X is 92.7 times larger

t(design.o) %*% design.o -> XpX
t(design.o) %*% y -> Xpy
chol(XpX) -> chol.o

determinant(chol.o)

b1 <- backsolve(chol.o,Xpy) # least squares solutions in two steps
b2 <- solve(XpX,Xpy)        # least squares estimates in one step
b3 <- backsolve(chol.o, forwardsolve(chol.o, Xpy),
                twice = FALSE) # in three steps
## checking that these three are indeed equal :
stopifnot(all.equal(b1, b2), all.equal(b2, b3))