SpatioTemporal (version 1.1.2)

makeCholBlock: Computations for Block Diagonal Matrices

Description

Provides block diagonal version of the base package functions chol, chol2inv, and backsolve. Computes the Cholesky factor, the matrix inverse and solves matrix equation systems for block diagonal matrices.

Usage

makeCholBlock(mat, n.blocks = 1,
    block.sizes = rep(dim(mat)[1]/n.blocks, n.blocks))
  invCholBlock(R, n.blocks = 1,
    block.sizes = rep(dim(R)[1]/n.blocks, n.blocks))
  solveTriBlock(R, B, n.blocks = 1,
    block.sizes = rep(dim(R)[1]/n.blocks, n.blocks),
    transpose = FALSE)

Arguments

mat

A block diagonal, square, positive definite matrix.

R

Upper right block diagonal Cholesky factor. The output from chol or makeCholBlock.

n.blocks

Number of diagonal blocks in mat (or R). Defaults to 1 (i.e. a full matrix) if neither n.blocks nor block.sizes given, o.w. it defaults to length(block.sizes)).

block.sizes

A vector of length n.blocks with the size of each of the diagonal blocks. If not given it will assume equal size blocks.

B

Vector or matrix containg the right hand side of the equations system to be solved; needs to be a multiple of dim(R)[1].

transpose

Transpose R before solving the equation system. Controls if we solve the equations system given by R*x = B or R'*x=B.

Value

• makeCholBlock gives the Cholesky factor and invCholBlock gives the inverse of the matrix mat. solveTriBlock gives to answer to the equation system.

Details

makeCholBlock computes the Cholesky factor of a block diagonal matrix using the block diagonal structure to speed up computations.

invCholBlock uses the Cholesky factor from makeCholBlock to compute the inverse of mat.

solveTriBlock solves equation systems based on the Cholesky factor, using the block diagonal structure to speed up computations (c.f. backsolve). The function solves equations of the form R*x = B, and R'*x = B with respect to x, where the transpose is controlled by the parameter transpose. Aplying the function twice solves mat*x=B, see the examples.

For all three functions the block diagonal structure of the matrix is defined by two input variables, the number of blocks n.blocks, and the size of each block block.sizes. The size of the matrices must match the total number of blocks, i.e. sum(block.sizes) must equal dim(mat).

The functions can be used for full matrices by setting the number of blocks to 1.

See Also

Other basic linear algebra: blockMult, crossDist, dotProd, norm2, sumLogDiag

Other block matrix functions: blockMult, calc.FX, calc.FXtF2, calc.iS.X, calc.mu.B, calc.tFX, calc.tFXF, calc.X.iS.X, makeSigmaB, makeSigmaNu

Examples

##create a matrix
mat <- cbind(c(1,0,0),c(0,2,1),c(0,1,2))
##define the number of blocks and block sizes
block.sizes <- c(1,2)
n.blocks <- length(block.sizes)

##Compute the Cholesky factor
R <- makeCholBlock(mat, n.blocks, block.sizes)
##and the matrix inverse
i.mat <- invCholBlock(R, n.blocks, block.sizes)
##compare to the alternative
i.mat-solve(mat)
if( max(abs(R-chol(mat))) > 1e-10 ){
    stop("makeCholBlock: Results not equal")
  }
  if( max(abs(i.mat-solve(mat))) > 1e-10 ){
    stop("invCholBlock: Results not equal")
  }
##define a B vector
B <- c(1,2,3)
##solve the equation system (we need n.x since B is not a matrix)
x1 <- solveTriBlock(R, B, n.blocks, block.sizes, tr=TRUE)
x2 <- solveTriBlock(R, x1, n.blocks, block.sizes, tr=FALSE)
print(x2)
##compare to the alternative
print(solve(mat,B))
range(x2-solve(mat,B))
if( max(abs(x2-solve(mat,B))) > 1e-10 ){
    stop("solveTriBlock: Results not equal")
  }
##compute the quadratic form B'*i.mat*B
norm2(x1)
##compare to the alternative
t(B) %*% i.mat %*% B
if( abs(norm2(x1) - t(B) %*% i.mat %*% B) > 1e-10 ){
    stop("solveTriBlock: Results not equal for norm")
  }