Provides block diagonal version of the base package functions
Computes the Cholesky factor, the matrix inverse and solves matrix equation systems for block diagonal matrices.
1 2 3 4 5 6 7 8
A block diagonal, square, positive definite matrix.
Number of diagonal blocks in
A vector of length
Upper right block diagonal Cholesky factor. The output from
Vector or matrix containg the right hand side of the equations
system to be solved; needs to be a multiple of
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
compute the inverse of
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
transpose. Aplying the function twice solves mat*x=B, see
For all three functions the block diagonal structure of the matrix is
defined by two input variables, the number of blocks
the size of each block
block.sizes. The size of the matrices must
match the total number of blocks, i.e.
The functions can be used for full matrices by setting the number of blocks to 1.
makeCholBlock gives the Cholesky factor and
invCholBlock gives the inverse of the matrix
solveTriBlock gives to answer to the equation system.
Johan Lindstrom and Adam Szpiro
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
##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) ##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)) ##compute the quadratic form B'*i.mat*B norm2(x1) ##compare to the alternative t(B) %*% i.mat %*% B
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.