mgcv (version 1.9-5)

bandchol: Cholesky decomposition of a band diagonal matrix

Description

Computes Cholesky decomposition of a (symmetric positive definite) band-diagonal matrix, A, or optionally a partial Cholesky factor if the matrix consists only of diagonals of a PD matrix without necessarily being PD. Solves or multiplies with the corresponding banded Cholesky factor R.

Usage

bandchol(B,partial=FALSE)
bandmult(B,x,trans=FALSE)
bandsolve(B,x,trans=FALSE)

Value

For bandchol, let R be the factor such that t(R)%*%R = A. R is upper triangular and if the rows of B contained the diagonals of A on entry, then what is returned is an n by k matrix containing the diagonals of R, packed as B was packed on entry. If B was square on entry, then R is returned directly. Returns NULL if matrix defined by B is not positive definite and partial=FALSE. See examples.

For bandmult or bandsolve a vector or matrix of the same dimension as x containing the required result.

Arguments

B

An n by k matrix containing the diagonals of the matrix A to be decomposed, or factor R. First row is leading diagonal, next is first super-diagonal, etc. Super-diagonals are zero padded at the end. Alternatively, for bandchol only, gives A directly, i.e. a square matrix with 2k-1 non zero diagonals (those from the lower triangle are not accessed).

partial

compute partial factor if B contains diagonals of a PD matrix, but the matrix with just these bands is not guaranteed PD. See details.

x

an n-vector or n row matrix to multiply by R (bandmult) or its inverse (bandsolve).

trans

if TRUE then transpose R before multiplication or solving.

Author

Simon N. Wood [email protected]

Details

For bandchol calls dpbtrf from LAPACK. The point of this is that it has \(O(k^2n)\) computational cost, rather than the \(O(n^3)\) required by dense matrix methods. bandmult calls BLAS routine dtbmv (repeatedly if x has several columns). bandsolve calls dtbtrs from LAPACK.

If partial=TRUE then bandchol performs a partial decomposition (calls compiled C code, not LAPACK). The diagonals in B should come from a PD matrix, but the matrix consisting of just those diagonals need not be PD. The leading diagonal must be positive. The factor computed is of a block diagonal approximation to the matrix whose bands are in B. The approximation's elements either match the target matrix exactly, or are set to zero. The construction runs the Cholesky decomposition until an imaginary diagonal factor element is encountered, at which point a new diagonal block is started, and all already computed elements beyond the previous block are set to zero. This is repeated until a full factor is produced. See examples.

References

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. and Sorensen, D., 1999. LAPACK Users' guide (Vol. 9). SIAM.

Examples

require(mgcv)
## simulate a banded diagonal matrix
n <- 7;set.seed(8)
A <- matrix(0,n,n)
sdiag(A) <- runif(n);sdiag(A,1) <- runif(n-1)
sdiag(A,2) <- runif(n-2)
A <- crossprod(A) 

## full matrix form...
bandchol(A)
R <- chol(A); R ## for comparison

## compact storage form...
B <- matrix(0,3,n)
B[1,] <- sdiag(A);B[2,1:(n-1)] <- sdiag(A,1)
B[3,1:(n-2)] <- sdiag(A,2)
RB <- bandchol(B); RB

## bandmult and bandsolve with a vector...
x <- runif(n)

drop(R%*%x) ## direct
bandmult(RB,x) ## banded
drop(t(R)%*%x)
bandmult(RB,x,trans=TRUE)

backsolve(R,x)
bandsolve(RB,x)
forwardsolve(t(R),x)
bandsolve(RB,x,trans=TRUE)

## bandmult and bandsolve with a matrix...
x <- matrix(runif(2*n),n,2)

drop(R%*%x)
bandmult(RB,x)
drop(t(R)%*%x)
bandmult(RB,x,trans=TRUE)

backsolve(R,x)
bandsolve(RB,x)
forwardsolve(t(R),x)
bandsolve(RB,x,trans=TRUE)

## Partial factor example
set.seed(0); n <- 7
B0 <- crossprod(matrix(runif(n*n),n,n)) ## PD matrix
b <- 3 ## super-diagonals to use
B <- matrix(0,b+1,n) ## packed storage for b+1 diags
B[1,] <- diag(B0)
for (i in 1:b) B[i+1,1:(n-i)] <- sdiag(B0,i)
bandchol(B) ## not PD
R <- bandchol(B,partial=TRUE) ## get partial factor
## pack back into full matrix for demo...
R0 <- diag(R[1,]); for (i in 1:b) sdiag(R0,i) <- R[i+1,1:(n-i)]
BR <- crossprod(R0)
## c.f. non zero crossprod(R0) and corresponding B0  
BR;B0*as.numeric(BR!=0)