块对角矩阵答案

【问题标题】：block diagonal matrix块对角矩阵
【发布时间】：2016-03-21 21:51:44
【问题描述】：

假设我有一个由许多矩阵组成的列表Z，我想从中构造一个块对角矩阵。例如：

[[1]]
             [,1]        [,2]         [,3]
[1,] 1.002500e+00 0.001930454 1.388794e-11
[2,] 1.930454e-03 1.002500000 1.930454e-03
[3,] 1.388794e-11 0.001930454 1.002500e+00

[[2]]
             [,1]        [,2]         [,3]
[1,] 1.002500e+00 0.001930454 1.388794e-11
[2,] 1.930454e-03 1.002500000 1.930454e-03
[3,] 1.388794e-11 0.001930454 1.002500e+00

我想创建一个块对角矩阵，我目前正在使用

block = bdiag(z)

但是，当列表中的矩阵数量很大时，bdiag 命令会很慢。从列表中构造块对角矩阵的快速简便方法是什么？

注意我的矩阵也是对称的，列表中的每个矩阵都有相似的维度。

【问题讨论】：

列表有多大？
@rawr 大约 200 个矩阵，每个 25*25
我尝试了 100,000 个 3x3 矩阵，不到一分钟。 200 个 25x25 矩阵的列表，例如，l <- rep(list(matrix(1:625, 25)), 200); x <- Matrix::bdiag(l) 几乎立即完成
@rawr 我的应用程序需要优化涉及 bdiag() 的函数，因此如果我的 bdiag 需要 1 秒，那么我的优化将花费大量时间。我正在尽可能地压缩代码。
它可能可以通过索引到稀疏矩阵来更快地做到这一点，但bdiag() 步骤似乎不太可能真的是限制因素......？

标签： r matrix symmetric

【解决方案1】：

我试图做一个更快的功能，bdiag2。

library(Matrix)

bdiag2 <- function(Ms){
  l <- length(Ms)
  N <- nrow(Ms[[1]])  
  i0 <- rep(1:N, times=N:1)
  s <- rep(seq(0,(l-1)*N,by=N),each=length(i0))
  i <- rep(i0,l) + s
  j0 <- 1:N -> j
  for(k in 1:N){
    j0 <- j0[-1]
    j <- c(j, j0)
  }
  j <- rep(j,l) + s
  idx <- t(upper.tri(Ms[[1]], diag = TRUE))
  x <- unlist(lapply(Ms, "[", idx))
  sparseMatrix(i, j, x = x, symmetric = TRUE)
}

# test & benchmark
N <- 5; l <- 200
Ms <- replicate(l, toeplitz(1:N), simplify = FALSE)

stopifnot(all(bdiag(Ms) == bdiag2(Ms)))

library(microbenchmark)
microbenchmark(
  bdiag = bdiag(Ms),
  bdiag2 = bdiag2(Ms)
)


> microbenchmark(
+   bdiag = bdiag(Ms),
+   bdiag2 = bdiag2(Ms)
+ )
Unit: milliseconds
   expr      min        lq      mean    median        uq       max neval cld
  bdiag 25.68078 27.393898 29.710474 28.566398 30.265242 54.393319   100   b
 bdiag2  1.34185  1.476838  1.693573  1.558724  1.769127  4.482054   100  a

【讨论】：