Ben Bolker
Ben Bolker

Reputation: 226192

construct new supermatrix from block matrices

How can I construct (in R) a matrix made of subcomponents that are matrices?

For example, starting from matrices

A <- matrix(1:9,nrow=3,ncol=3)
B <- matrix(5:10,nrow=2,ncol=3)
C <- matrix(11:20,nrow=2,ncol=5)

I want to construct a block matrix

A 0
B C

where 0 represents a zero-filled block with the appropriate dimensions.

There are other questions on SO about constructing block-diagonal matrices (Matrix::bdiag is very good for this), but I can't find one that answers this question.

(I'm posting this question because I had just about finished answering it when it was deleted by its original poster ...)

Upvotes: 1

Views: 223

Answers (2)

flodel
flodel

Reputation: 89057

I tried writing a general purpose function. The usage is similar to matrix() but the first argument is a list of matrices (or vectors that will be recycled). It does not have all the bells and whistles (dimnames, byrow) but it is a decent start. I wouldn't be surprised to find out a better and more complete function already exists in a package but at least it was a fun exercise:

supermatrix <- function(list.of.mat, nrow = 1L, ncol = 1L) {
   stopifnot(length(list.of.mat) == nrow * ncol)
   is.mat <- vapply(list.of.mat, is.matrix, logical(1L))
   is.vec <- vapply(list.of.mat, is.vector, logical(1L))
   if (any(!is.mat & !is.vec)) stop("the list items must be matrices or vectors")
   is.mat.mat <- matrix(is.mat, nrow, ncol)
   if (any(rowSums(is.mat.mat) == 0L))
      stop("we need at least one matrix per super row")
   if (any(colSums(is.mat.mat) == 0L))
      stop("we need at least one matrix per super column")
   na.mat <- matrix(NA, nrow, ncol)
   nrow.mat <- replace(na.mat, is.mat, vapply(list.of.mat[is.mat], nrow, integer(1L)))
   ncol.mat <- replace(na.mat, is.mat, vapply(list.of.mat[is.mat], ncol, integer(1L)))
   is.not.uniq <- function(x) length(table(x)) > 1L
   if (any(apply(nrow.mat, 1, is.not.uniq))) stop("row dim mismatch")
   if (any(apply(ncol.mat, 2, is.not.uniq))) stop("col dim mismatch")
   nrow.vec <- rowMeans(nrow.mat, na.rm = TRUE)
   ncol.vec <- colMeans(ncol.mat, na.rm = TRUE)
   nrow.mat <- matrix(nrow.vec, nrow, ncol, byrow = FALSE)
   ncol.mat <- matrix(ncol.vec, nrow, ncol, byrow = TRUE)
   all.mat <- Map(matrix, list.of.mat, nrow.mat, ncol.mat)
   i1.idx <- unlist(Map(rep, row(na.mat), lapply(all.mat, length)))
   j1.idx <- unlist(Map(rep, col(na.mat), lapply(all.mat, length)))
   i2.idx <- unlist(lapply(all.mat, row))
   j2.idx <- unlist(lapply(all.mat, col)) 
   o.idx <- order(j1.idx, j2.idx, i1.idx, i2.idx)
   matrix(unlist(all.mat)[o.idx], sum(nrow.vec), sum(ncol.vec))
}

Example usage:

A <- matrix(1:9,nrow=3,ncol=3)
B <- matrix(5:10,nrow=2,ncol=3)
C <- matrix(11:20,nrow=2,ncol=5)

supermatrix(list(A, B, 0, C), 2, 2)
supermatrix(list(A, B, A, 1, 0, C, 2, C), 4, 2)

Upvotes: 4

Ben Bolker
Ben Bolker

Reputation: 226192

We need a zero matrix that will have compatible dimensions with A and C:

z <- matrix(0,nrow=nrow(A),ncol=ncol(C))

Now we just use rbind() and cbind():

rbind(cbind(A,z),cbind(B,C))

Upvotes: 0

Related Questions