Block-diagonal binding of matrices
adiag from a package magic does what you want:
library(magic)
adiag(a,b)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 3 5 0 0
[2,] 2 4 6 0 0
[3,] 0 0 0 7 9
[4,] 0 0 0 8 10
Alternatively, you could use a package Matrix and function bdiag
library(Matrix)
bdiag(a,b)
4 x 5 sparse Matrix of class "dgCMatrix"
[1,] 1 3 5 . .
[2,] 2 4 6 . .
[3,] . . . 7 9
[4,] . . . 8 10
that returns a sparse matrix and which might be more efficient. Use as.matrix(bdiag(a,b)) to get a regular one.
Build a block diagonal matrix from a known matrix
We can use bdiag
library(Matrix)
bdiag(replicate(3, B, simplify = FALSE))
#6 x 6 sparse Matrix of class "dgCMatrix"
#[1,] 1 4 . . . .
#[2,] 3 5 . . . .
#[3,] . . 1 4 . .
#[4,] . . 3 5 . .
#[5,] . . . . 1 4
#[6,] . . . . 3 5
Can we wrapped in a function
fdiag <- function(mat, n) {
bdiag(replicate(n, mat, simplify = FALSE))
}
fdiag(B, 3)
How to create a diagonal matrix with blocks in the diagonal with each block a different matrix?
We can use bdiag from Matrix after placing the matrices in a list to get a sparseMatrix.
library(Matrix)
sM <- bdiag(mget(paste0("A", 1:5)))
sM
# 10 x 10 sparse Matrix of class "dgCMatrix"
# [1,] 1 3 . . . . . . . .
# [2,] 2 4 . . . . . . . .
# [3,] . . 5 7 . . . . . .
# [4,] . . 6 8 . . . . . .
# [5,] . . . . 9 11 . . . .
# [6,] . . . . 10 12 . . . .
# [7,] . . . . . . 13 15 . .
# [8,] . . . . . . 14 16 . .
# [9,] . . . . . . . . 17 19
#[10,] . . . . . . . . 18 20
If needed, the sparseMatrix can be converted to a matrix
as.matrix(sM)
data
A1 <- matrix(1:4, 2, 2)
A2 <- matrix(5:8, 2, 2)
A3 <- matrix(9:12, 2, 2)
A4 <- matrix(13:16, 2, 2)
A5 <- matrix(17:20, 2, 2)
How could I build a function that extracts the diagonal block matrices of a larger one in R
I think a clearer solution is to reset the dimensions and then let R do the index calculations for you:
err_cov <- function(x){
m <- nrow(x)
n <- ncol(x)
#compute the full error covariance matrix as the inner product
#of vec(x) and its transpose
vec <- as.vector(x)
omega <- tcrossprod(vec)
dim(omega) <- c(n,m,n,m)
sigmas <- list()
for (j in 1:m)
sigmas[[j]] <- omega[,j,,j]
return(sigmas)
}
Here is an example:
> x
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> tcrossprod(vec)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 2 3 4 5 6
[2,] 2 4 6 8 10 12
[3,] 3 6 9 12 15 18
[4,] 4 8 12 16 20 24
[5,] 5 10 15 20 25 30
[6,] 6 12 18 24 30 36
> err_cov(x)
[[1]]
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 4 6
[3,] 3 6 9
[[2]]
[,1] [,2] [,3]
[1,] 16 20 24
[2,] 20 25 30
[3,] 24 30 36
Repeating a block matrix many times in diagonal part of a block matrix, with the off-diagonal blocks all zero matrices?
1) kronecker If M is your matrix and k is the number of times you want it repeated then:
kronecker(diag(k), M)
For example,
M <- matrix(0.5, 5, 5) + diag(0.5, 5)
k <- 2
kronecker(diag(k), M)
giving:
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1.0 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.0 0.0
[2,] 0.5 1.0 0.5 0.5 0.5 0.0 0.0 0.0 0.0 0.0
[3,] 0.5 0.5 1.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
[4,] 0.5 0.5 0.5 1.0 0.5 0.0 0.0 0.0 0.0 0.0
[5,] 0.5 0.5 0.5 0.5 1.0 0.0 0.0 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 1.0 0.5 0.5 0.5 0.5
[7,] 0.0 0.0 0.0 0.0 0.0 0.5 1.0 0.5 0.5 0.5
[8,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 1.0 0.5 0.5
[9,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.5 1.0 0.5
[10,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.5 0.5 1.0
1a) %x% The last line of code can alternately be written as:
diag(k) %x% M
2) Matrix::bdiag Another possibility if you want to save space is to create a sparse matrix of class "dgMCatrix". It does not store the zero values. See ?bdiag :
library(Matrix)
bdiag(replicate(k, M, simplify = FALSE))
giving:
10 x 10 sparse Matrix of class "dgCMatrix"
[1,] 1.0 0.5 0.5 0.5 0.5 . . . . .
[2,] 0.5 1.0 0.5 0.5 0.5 . . . . .
[3,] 0.5 0.5 1.0 0.5 0.5 . . . . .
[4,] 0.5 0.5 0.5 1.0 0.5 . . . . .
[5,] 0.5 0.5 0.5 0.5 1.0 . . . . .
[6,] . . . . . 1.0 0.5 0.5 0.5 0.5
[7,] . . . . . 0.5 1.0 0.5 0.5 0.5
[8,] . . . . . 0.5 0.5 1.0 0.5 0.5
[9,] . . . . . 0.5 0.5 0.5 1.0 0.5
[10,] . . . . . 0.5 0.5 0.5 0.5 1.0
2b) Diagonal or to create a sparse matrix of class "dgTMatrix" :
Diagonal(k) %x% M
Combine matrices at corner in R?
We could use bdiag
library(Matrix)
as.matrix(bdiag(m1, m2))
-output
[,1] [,2] [,3] [,4]
[1,] 1.0 0.5 0.00 0.00
[2,] 0.5 1.0 0.00 0.00
[3,] 0.0 0.0 0.50 0.25
[4,] 0.0 0.0 0.25 0.50
construct new supermatrix from block matrices
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)
vector into block matrix
Instead of manually doing the split, we can use %/%
k <- 3
lst <- split(a, (seq_along(a)-1)%/%k + 1)
do.call(rbind, lapply(split(lst, (seq_along(lst)-1) %/% k + 1),
function(x) do.call(cbind, lapply(x, function(y) diag(y)))))
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
# [1,] 1 0 0 2 0 0 3 0 0
# [2,] 0 1 0 0 2 0 0 3 0
# [3,] 0 0 1 0 0 2 0 0 3
# [4,] 2 0 0 4 0 0 6 0 0
# [5,] 0 2 0 0 4 0 0 6 0
# [6,] 0 0 2 0 0 4 0 0 6
# [7,] 3 0 0 6 0 0 9 0 0
# [8,] 0 3 0 0 6 0 0 9 0
# [9,] 0 0 3 0 0 6 0 0 9
generating efficient sub block matrices to avoid duplicated code
Given the periodic nature of the data the submatrices and the whole block matrix can be filled via the PAD= argumant to the RESHAPE intrinsic.
program blox
implicit none
integer m
integer n
integer, allocatable :: zero(:), sp3(:,:,:)
integer, allocatable :: b1(:,:), b2(:,:), b3(:,:)
integer, allocatable :: iKE(:,:)
character(80) fmt
m = 4
n = m**2
zero = reshape([integer::],[m+1],pad=[0])
b1 = reshape([integer::],[m,m],pad=[-4,1,zero(1:m-2),1])
b2 = reshape([integer::],[m,m],pad=[1,zero(1:m)])
b3 = reshape([integer::],[m,m],pad=zero(1:m+1))
sp3 = reshape([integer::],[m,m,m-2],pad=b3)
iKE = reshape(reshape([integer::],[m,m,m,m],pad=[b1,b2,sp3,b2],order=[1,3,2,4]),[n,n])
write(fmt,'(*(g0))') '(',m,'(',m,'(',m,'i3:1x)/))'
write(*,fmt) transpose(iKE)
end program blox
Notice how strings of zeros are created by padding out an empty array with zeros (zero) and then arrays with periodic data are filled by padding out an empty array array with a single period of data (bl, b2, and b3). Then a matrix consisting of copies of a block matrix is created by padding out an empty array with the block (sp3). Finally a block-periodic matrix is created by padding out an empty array with the sequence of blocks. The resulting matrix has to be read out in cross-order and then reshaped into the right dimensions (iKE).
Output:
-4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0
1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0 0
0 1 -4 1 0 0 1 0 0 0 0 0 0 0 0 0
0 0 1 -4 0 0 0 1 0 0 0 0 0 0 0 0
1 0 0 0 -4 1 0 0 1 0 0 0 0 0 0 0
0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0 0
0 0 1 0 0 1 -4 1 0 0 1 0 0 0 0 0
0 0 0 1 0 0 1 -4 0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0 -4 1 0 0 1 0 0 0
0 0 0 0 0 1 0 0 1 -4 1 0 0 1 0 0
0 0 0 0 0 0 1 0 0 1 -4 1 0 0 1 0
0 0 0 0 0 0 0 1 0 0 1 -4 0 0 0 1
0 0 0 0 0 0 0 0 1 0 0 0 -4 1 0 0
0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1 0
0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4 1
0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 -4
Surprising that you really need the explicit form of this matrix, though. Usually you see this kind of object handled more indirectly via sparse matrix techniques.
Related Topics
Data.Frame Without Ruining Column Names
Plotting with Ggplot2: "Error: Discrete Value Supplied to Continuous Scale" on Categorical Y-Axis
Create End of the Month Date from a Date Variable
How to Select Last N Observation from Each Group in Dplyr Dataframe
Forward and Backward Fill Data Frame in R
Code Chunk Font Size in Rmarkdown with Knitr and Latex
Reading Rdata File with Different Encoding
Using Different Scales as Fill Based on Factor
Finding Out Which Functions Are Called Within a Given Function
Count Number of Rows Matching a Criteria
What Is the Meaning of the Dollar Sign "$" in R Function()
Convert from Billion to Million and Vice Versa