高效地构建由乘积子矩阵组成的矩阵

Question

我希望构建一个矩阵，如下所示（其中A是实数矩阵，I是单位矩阵）：

我不想使用for循环。我尝试了以下方法：

sequence = 1:T
sapply(sequence, function(i) matrix(A%^%(i-1)))

但是失败了。我考虑在图片中创建第一列矩阵，然后沿着对角线继续复制计算出的矩阵，但我不知道如何实现。

编辑: 对我造成的麻烦非常抱歉。这是我在一个快速而简单的for循环中寻找的内容。

library("expm")

n<-5

A<-matrix(1, 2, 2)
output <- matrix(0, 5*2, 5*2)

for (i in 1:5) {
  for (j in i:1) {
    output[(2*(i-1)+1):(2*i),(2*(j-1)+1):(2*j)] = A %^% (i-j)
  }
}

- user241879

---

快速而简单的循环通常是最好的选择... - user20650

PS：您应该将解决方案发布在答案部分，而不是作为对问题的编辑。 - user20650

问题是，我想在一个时间敏感的模拟研究中运行上面的循环 - 所以我必须尽可能快地编码。上面的解决方案只是为了举例说明我想做什么。一方面，存在这样的缺点，即对于任意自然数i，A^i的计算会一遍又一遍地进行。我在最近的答案中纠正了这个问题。但我认为使用类似apply的函数可以使过程更快。 - user241879

Answer 1

我刚想到另外一件事情：

library("expm")

n<-5
A<-matrix(1, 2, 2)

output <- matrix(0, 5*2, 5*2)
output2 <- matrix(0, 5*2, 5*2)

start.time <- Sys.time()
for(dummy in 1:2000){
  for (i in 1:5) {
    for (j in i:1) {
      output[(2*(i-1)+1):(2*i),(2*(j-1)+1):(2*j)] = A %^% (i-j)
    }
  }
}
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken

start.time <- Sys.time()


for(dummy in 1:2000){
  for (i in 1:5) {
    output2[(2*(i-1)+1):(2*i),1:2] = A %^% (i-1)
  }

  rowLength = dim(output2)[1]

  for (i in 2:5) {
    output2[(2*i-1):rowLength,(2*i-1):(2*i)] = output2[1:(rowLength-(2*(i-1))),1:2]
  }
}
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken

这会稍微快一点，但我认为通过一些内置函数的魔法，结果可以变得更好。也许有人能帮忙？

- user241879

Answer 2

你是否看过文档https://stat.ethz.ch/R-manual/R-devel/library/base/html/lower.tri.html中介绍的lower.tri(x, diag = TRUE)函数？

我已经制作了一个示例，你可以在这里应用它。

具体步骤如下：

A <- matrix(1:25, 5, 5)  # construct the initial matrix
B <- A * lower.tri(A)    # convert A to a lower triangular matrix of booleans (with TRUE below diagonal, FALSE from diagonal and up)
B <- B * A + diag(5)     # construct the final result
B

- Kent Munthe Caspersen

---

谢谢你的回答。我刚看了一下，但是我不知道如何使用它来设置我的矩阵的行相关指数。另外，我尝试使用sapply，但是我得到了一个混乱的向量。 - user241879

@user241879，我刚刚附上了一个应用该函数的示例。 - Kent Munthe Caspersen

@user20650 我提供的示例使用正方形矩阵。你能举个例子说明输出结果与预期不同吗？ - Kent Munthe Caspersen

抱歉，我删除了评论...我的初步想法是A必须具有以下数量的元素才能填充lower.tri... n_lt <- function(n) 1/2 * n*(n+1) ; n_lt(1:5)（即如果标识为5维，则lower.tri有15个元素，因此A必须有15个元素，这显然对于方阵不可能发生）。但是，我认为我误解了问题...它似乎在对角线上重复了完整的矩阵A，并且A的幂等等，而不仅仅是元素。 - user20650

没问题，我错过了指数，所以我的答案不能解决问题。我把它留在这里，以防一些思路可以借鉴。 - Kent Munthe Caspersen

谢谢您的时间。我尝试使用lower.trig函数进行操作，但无法弄清楚如何使用它来获得我想要的最终结果。我还更新了我的问题以展示我确切的需求。无论如何，还是非常感谢您！ - user241879