在R中创建块循环矩阵



我试图在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

替代方法

这是一种使用kroneckerReduce 的高效方法

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

相关内容

  • 没有找到相关文章

最新更新