我试图在R中创建一个块循环矩阵。下面给出了块循环矩阵的结构。
C0 C1 ... Cn-1
Cn-1 C0 C1 ... Cn-2
Cn-2 Cn-1 .... Cn-3
and so on
我有块
C0 .... Cn-1
创建矩阵的最简单方法是什么。是否有可用的功能?
感谢您提出一个具有挑战性的问题!这是一个求和矩阵与子对角线和超对角线的kronecker乘积的解决方案。
样本数据,矩阵列表:
C <- lapply(1:3, matrix, nrow = 2, ncol = 2)
我的解决方案:
bcm <- function(C) {
require(Matrix)
n <- length(C)
Reduce(`+`, lapply((-n+1):(n-1),
function(i) kronecker(as.matrix(bandSparse(n, n, -i)),
C[[1 + (i %% n)]])))
}
bcm(C)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 1 1 3 3 2 2
# [2,] 1 1 3 3 2 2
# [3,] 2 2 1 1 3 3
# [4,] 2 2 1 1 3 3
# [5,] 3 3 2 2 1 1
# [6,] 3 3 2 2 1 1
我不知道这是否特别有效,但当我解释你的问题时,它会满足你的需求。
rotList <- function(L,n) {
if (n==0) return(L)
c(tail(L,n),head(L,-n))
}
rowFun <- function(n,matList) do.call(rbind,rotList(matList,n))
bcMat <- function(matList) {
n <- length(matList)
do.call(cbind,lapply(0:(n-1),rowFun,matList))
}
示例:
bcMat(list(diag(3),matrix(1:9,nrow=3),matrix(4,nrow=3,ncol=3)))
我想您要查找的是lgcp
包中的circulant.matrix
。
如果x是一个矩阵,其列是块循环矩阵,则此函数返回块循环感兴趣的矩阵。
例如
x <- matrix(1:8,ncol=4)
circulant(x)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,] 1 2 3 4 5 6 7 8
# [2,] 2 1 4 3 6 5 8 7
# [3,] 7 8 1 2 3 4 5 6
# [4,] 8 7 2 1 4 3 6 5
# [5,] 5 6 7 8 1 2 3 4
# [6,] 6 5 8 7 2 1 4 3
# [7,] 3 4 5 6 7 8 1 2
# [8,] 4 3 6 5 8 7 2 1
替代方法
这是一种使用kronecker
和Reduce
的高效方法
bcirc <- function(list.blocks){
P <- lapply(seq_along(list.blocks), function(x,y) x ==y, x = circulant(seq_along(list.blocks)))
Reduce('+',Map(P = P, A=list.blocks, f = function(P,A) kronecker(P,A)))
}
与@flodel和@Ben Bolker 进行基准测试
lbirary(microbenchmark)
microbenchmark(bcm(C), bcirc(C), bcMat(C))
Unit: microseconds
expr min lq median uq max neval
bcm(C) 10836.719 10925.7845 10992.8450 11141.1240 21622.927 100
bcirc(C) 444.983 455.7275 479.5790 487.0370 569.105 100
bcMat(C) 288.558 296.4350 309.8945 348.4215 2190.231 100
你想要这样的东西吗?
> vec <- 1:4
> sapply(rev(seq_along(vec)),function(x) c(tail(vec,x),head(vec,-x)) )
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 2 3 4 1
[3,] 3 4 1 2
[4,] 4 1 2 3