Linear Equation Solving for Sparse Matrices

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.

# S4 method for spam
solve(a, b, Rstruct=NULL, ...)
# S4 method for spam
backsolve(r, x, ...)
# S4 method for spam
forwardsolve(l, x, ...)
# S4 method for spam
chol2inv(x, ...)

Arguments

a

symmetric positive definite matrix of class spam or a Cholesky factor as the result of a chol call.

l,r

object of class spam or 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.

Rstruct

the Cholesky structure of a.

...

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.

In case a is from a chol call, then solve is an efficient way to calculate backsolve(a, forwardsolve( t(a), b)).

However, for a.spam and a.mat from a chol call with a sparse and ordinary matrix, note that forwardsolve( a.mat, b, transpose=T, upper.tri=T) is equivalent to forwardsolve( t(a.mat), b) and backsolve(a.spam, forwardsolve(a.spam, b, transpose=T, upper.tri=T)) yields the desired result. But backsolve(a.spam,forwardsolve(t(a.spam), resid)) is wrong because t(a.spam) is a spam and not a spam.chol.NgPeyton object.

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.

Notice that it is more efficient to solve successively the linear equations (both triangular solves) than to implement these in the Fortran code.

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.

chol2inv(x) is a faster way to solve(x).

Note

There is intentionally no S3 distinction between the classes spam and spam.chol.method.

References

See references in chol.

See also

chol.spam and ordering.

Author

Reinhard Furrer, based on Ng and Peyton (1993) Fortran routines

Examples

# Generate multivariate form a covariance inverse:
# (usefull for GRMF)
set.seed(13)
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)
#> [1] 1.049394

# compare with:
mvsample <- backsolve( chol(as.matrix( Sigmainv)), iidsample, n)
   #### ,n as patch
norm( var(t(mvsample)) - Sigma)
#> [1] 1.040029



# 'solve' step by step:
b <- rnorm( n)
R <- chol(Sigmainv)
norm( backsolve( R, forwardsolve( R, b))-
      solve( Sigmainv, b) )
#> [1] 0
norm( backsolve( R, forwardsolve( R, diag(n)))- Sigma )
#> [1] 0


# 'update':
R1 <- update( R, Sigmainv + diag.spam( n))