bdiag {Matrix} | R Documentation |

Build a block diagonal matrix given several building block matrices.

bdiag(...) .bdiag(lst)

`...` |
individual matrices or a |

`lst` |
non-empty |

For non-trivial argument list, `bdiag()`

calls `.bdiag()`

.
The latter maybe useful to programmers.

A *sparse* matrix obtained by combining the arguments into a
block diagonal matrix.

The value of `bdiag()`

inheris from class
`CsparseMatrix`

, whereas
`.bdiag()`

returns a `TsparseMatrix`

.

This function has been written and is efficient for the case of relatively few block matrices which are typically sparse themselves.

It is currently *inefficient* for the case of many small dense
block matrices.
For the case of *many* dense *k * k* matrices,
the `bdiag_m()`

function in the ‘Examples’ is an order of
magnitude faster.

Martin Maechler, built on a version posted by Berton Gunter to
R-help; earlier versions have been posted by other authors, notably
Scott Chasalow to S-news. Doug Bates's faster implementation builds
on `TsparseMatrix`

objects.

`Diagonal`

for constructing matrices of
class `diagonalMatrix`

, or `kronecker`

which also works for `"Matrix"`

inheriting matrices.

`bandSparse`

constructs a *banded* sparse matrix from
its non-zero sub-/super - diagonals.

Note that other CRAN **R** packages have own versions of `bdiag()`

which return traditional matrices.

bdiag(matrix(1:4, 2), diag(3)) ## combine "Matrix" class and traditional matrices: bdiag(Diagonal(2), matrix(1:3, 3,4), diag(3:2)) mlist <- list(1, 2:3, diag(x=5:3), 27, cbind(1,3:6), 100:101) bdiag(mlist) stopifnot(identical(bdiag(mlist), bdiag(lapply(mlist, as.matrix)))) ml <- c(as(matrix((1:24)%% 11 == 0, 6,4),"nMatrix"), rep(list(Diagonal(2, x=TRUE)), 3)) mln <- c(ml, Diagonal(x = 1:3)) stopifnot(is(bdiag(ml), "lsparseMatrix"), is(bdiag(mln),"dsparseMatrix") ) ## random (diagonal-)block-triangular matrices: rblockTri <- function(nb, max.ni, lambda = 3) { .bdiag(replicate(nb, { n <- sample.int(max.ni, 1) tril(Matrix(rpois(n*n, lambda=lambda), n,n)) })) } (T4 <- rblockTri(4, 10, lambda = 1)) image(T1 <- rblockTri(12, 20)) ##' Fast version of Matrix :: .bdiag() -- for the case of *many* (k x k) matrices: ##' @param lmat list(<mat1>, <mat2>, ....., <mat_N>) where each mat_j is a k x k 'matrix' ##' @return a sparse (N*k x N*k) matrix of class \code{"\linkS4class{dgCMatrix}"}. bdiag_m <- function(lmat) { ## Copyright (C) 2016 Martin Maechler, ETH Zurich if(!length(lmat)) return(new("dgCMatrix")) stopifnot(is.list(lmat), is.matrix(lmat[[1]]), (k <- (d <- dim(lmat[[1]]))[1]) == d[2], # k x k all(vapply(lmat, dim, integer(2)) == k)) # all of them N <- length(lmat) if(N * k > .Machine$integer.max) stop("resulting matrix too large; would be M x M, with M=", N*k) M <- as.integer(N * k) ## result: an M x M matrix new("dgCMatrix", Dim = c(M,M), ## 'i :' maybe there's a faster way (w/o matrix indexing), but elegant? i = as.vector(matrix(0L:(M-1L), nrow=k)[, rep(seq_len(N), each=k)]), p = k * 0L:M, x = as.double(unlist(lmat, recursive=FALSE, use.names=FALSE))) } l12 <- replicate(12, matrix(rpois(16, lambda = 6.4), 4,4), simplify=FALSE) dim(T12 <- bdiag_m(l12))# 48 x 48 T12[1:20, 1:20]

[Package *Matrix* version 1.2-17 Index]