5

我正在尝试在 R 中创建一个块循环矩阵。块循环矩阵的结构如下所示。

C0 C1 ... Cn-1
Cn-1 C0 C1 ... Cn-2
Cn-2 Cn-1 .... Cn-3

and so on

我有积木

C0 .... Cn-1

创建矩阵的最简单方法是什么。有没有已经可用的功能？

4

4 回答 4

6

感谢您提出具有挑战性的问题！这是一个将矩阵的克罗内克乘积与子对角线和超对角线相加的解决方案。

样本数据，矩阵列表：

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
于 2013-07-04T04:06:28.243 回答
4

我不知道这是否特别有效，但是当我解释您的问题时，它可以满足您的要求。

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)))
于 2013-07-04T03:56:48.110 回答
3

我认为您正在寻找的是circulant.matrix来自lgcp包装的东西。

如果 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

替代方法

这是一种非常低效的方法，使用kroneckerandReduce

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
于 2013-07-04T03:55:59.390 回答
0

你正在寻找这样的东西吗？

> 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
于 2013-07-04T03:25:25.480 回答