spam (version 0.11)

spam solve: Linear Equation Solving for Sparse Matrices

Description

backsolve and forwardsolve solve a system of linear equations where the coefficient matrix is upper or lower triangular. solve solves a linear system or computes the inverse of a matrix if the right-hand-side is missing.

Usage

solve(a, b, ...)
backsolve(r, x, ...)forwardsolve(l, x, ...)

Arguments

a

symmetric positive definite matrix of class spam.

l,r

object of class spam.chol.method returned by the function chol.

x,b

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

...

further arguments passed to or from other methods, see Details below.

Details

We can solve A %*% x = b by first computing the Cholesky decomposition A = t(R)%*%R), then solving t(R)%*%y = b for y, and finally solving R%*%x = y for x. solve combines chol, a Cholesky decomposition of a symmetric positive definite sparse matrix, with forwardsolve and then backsolve.

forwardsolve and backsolve solve a system of linear equations where the coefficient matrix is lower (forwardsolve) or upper (backsolve) triangular. Usually, the triangular matrix is result from a chol call and it is not required to transpose it for forwardsolve. Note that arguments of the default methods k, upper.tri and transpose do not have any effects here.

If the right-hand-side in solve is missing it will compute the inverse of a matrix. For details about the specific Cholsesky decomposition, see chol.

Recall that the Cholesky factors are from ordered matrices.

References

See also references in chol. Ng, E. G. and B. W. Peyton (1993), "Block sparse Cholesky algorithms on advanced uniprocessor computers", SIAM J. Sci. Comput., 14, pp. 1034-1056.

See Also

det, chol and ordering.

Examples

# Generate multivariate form a covariance inverse:
# (usefull for GRMF)
n <- 25    # dimension
N <- 1000  # sample size
Sigmainv <- .25^abs(outer(1:n,1:n,"-"))
Sigmainv <- as.spam( Sigmainv, eps=1e-4)


Sigma <- solve( Sigmainv)  # for verification 
iidsample <- array(rnorm(N*n),c(n,N))

mvsample <- backsolve( chol(Sigmainv), iidsample)
norm( var(t(mvsample)) - Sigma, type="HS")

# compare with:
mvsample <- backsolve( chol(as.matrix( Sigmainv)), iidsample)
norm( var(t(mvsample)) - Sigma, type="HS")



# 'solve' step by step:
b <- rnorm( n)
R <- chol(Sigmainv)
norm( backsolve( R, forwardsolve( R, b))-
      solve( Sigmainv, b),type="HS") 
norm( backsolve( R, forwardsolve( R, diag(n)))- Sigma,type="HS")