这个问题可能会引起统计学家和数据科学家的兴趣，因为它提出了多维数组处理问题。

循环不是问题。 您的目标可能是高效地执行操作，不管它是如何完成的，并且可能以一种明显可并行的方式来完成。

利用 R 管理多维数组的能力。 （它从 FORTRAN 继承了这一点，即使在其原始设计中也可以处理多达七个维度。）将数据组整体移动到一个数组中，然后根据需要对其进行整形。即使需要进行一些循环，这些本机操作也会很快。

尽管它使用短的高级循环（每个循环执行次），但以下解决方案将创建一个和个块的“block Toeplitz”数组——一个包含四百万个元素的数组——在一个半秒和少量使用开销存储。它基于观察到可以从字符串构造 Toeplitz 模式$n+1$$n=999$$2\times 2$

通过取最后元素形成输出的第一列，将这些元素向左移动形成第二列输出，依此类推，直到输出的最后一列由前元素形成这个字符串。$n+1$$n+1$

该代码用于n表示（因为索引从 1 而不是 0 开始）。矩阵串本身，作为一个数组，存储在一个临时数组中。$n+1$R$k\times (2n-1)$strip

最后，通过“爆破”输出列的条带部分创建数组被置换以使其值以正确的顺序存储（我懒得一开始就制定出最佳存储配置），然后重铸为数组（不涉及数据移动）。这些按列操作的每个操作都是高度矢量化的，并且能够进行大规模并行。$k\times k\times n\times n$X$kn\times kn$

的初始值U包含用于检查输出是否正确的明显模式。

#
# Create n square matrices of dimension k by k.
#
U <- list(matrix(1:4, 2), matrix(5:8, 2), matrix(9:12, 2))
#U <- lapply(1:1000, function(i) matrix(-3:0 + 4*i, 2))

system.time({
  k <- min(unlist(lapply(U, dim)))
  n <- length(U)
  #
  # Create the "strip".
  #
  strip <- array(NA, dim=c(k,k,2*n-1))
  for (i in 1:n) strip[,,i] <- U[[n+1-i]]
  if (n > 1) for (i in 2:n) strip[,,n+i-1] <- t(U[[i]])
  #
  # Assemble into "block-Toeplitz" form.
  #
  X <- array(NA, dim=c(k,k,n,n))
  #
  # Blast the strip across X.
  #
  for (i in 1:n) X[,,,i] <- strip[,,(n+1-i):(2*n-i)]
  X <- matrix(aperm(X, c(1,3,2,4)), n*k)
})

以下函数将块列表作为参数。但是，这可能不是最好的解决方案：仍然存在一个循环。而且它需要更多的工作，因为它不做块下对角线的转置（在我的例子中，我有对称矩阵）。

toeplitz.block <- function(blocks) {
    l <- length(blocks)
    m.str <- toeplitz(1:l)

    res <- lapply(1:l,function(k) {
        res <- matrix(0,ncol=ncol(m.str),nrow=nrow(m.str))
        res[m.str == k] <- 1
        res %x% blocks[[k]]
    })

    Reduce("+",res)
}