Linear Equation Solving for Sparse Matrices
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 and will
compute the inverse of a matrix if the right-hand-side is missing.
chol(x, ...)
## S4 method for signature 'matrix.csr.chol'
backsolve(r, x, k = NULL, upper.tri = NULL,
transpose = NULL, twice = TRUE, ...)
forwardsolve(l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)
solve(a, b, ...)a |
symmetric positive definite matrix of class |
r |
object of class |
l |
object of class |
x,b |
vector(regular matrix) of right-hand-side(s) of a system of linear equations. |
k |
inherited from the generic; not used here. |
upper.tri |
inherited from the generic; not used here. |
transpose |
inherited from the generic; not used here. |
twice |
Logical flag: If true backsolve solves twice, see below. |
... |
further arguments passed to or from other methods. |
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.
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.
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. Sci. Comput., 14, pp. 1034-1056.
slm for sparse version of lm
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
t(design.o) %*% design.o -> XpX
t(design.o) %*% y -> Xpy
chol(XpX) -> 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))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.