在 R 中填充矩阵的更快方法答案

【问题标题】：Faster way of filling a matrix in R在 R 中填充矩阵的更快方法
【发布时间】：2021-07-19 05:08:59
【问题描述】：

我想在 R 中填充一个矩阵，但每一列都必须有一个迭代向下移动的向量。所以从某种意义上说，它将是一个下三角矩阵。我的努力是这样的：


x = c(3,4,8,9)
E <- matrix(0,length(x),length(x));E
for (i in 1:nrow(E)){
  E[i,1]=x[i]
}
E
for (i in 2:nrow(E)){
  for (j in 2:ncol(E)) {
    E[i,2] =x[i-1] }  }
E
for (i in 3:nrow(E)){
  for (j in 3:ncol(E)) {
    E[i,3] =x[i-2] }  }
E
for (i in 4:nrow(E)){
  for (j in 4:ncol(E)) {
    E[i,4] =x[i-3] }  }
E

每次从向量中删除一个元素。但是有没有一种更快的方法可以用更少的 for 循环和向量的 n 长度而不是 4 来完成它，作为一般化？

【问题讨论】：

什么是3 只存在于一个单元格中？这似乎违反直觉......

标签： r loops for-loop matrix

【解决方案1】：

对不起，我无法抗拒。这是另一种基本方法：

x <- c(3,4,8,9)
n <- length(x)
E <- diag(rep(x[1], n))
j <- unlist(sapply(length(x):2, function(i) x[2:i]))
E[lower.tri(E)] <- j

添加到 Rui 的基准代码中，我们得到：

【讨论】：

不错的基准测试。投票您的贡献。我还在回答中添加了一个选项。

【解决方案2】：

我认为如果您将此代码添加到基准测试中会很有趣

TIC <- function(x) {
  E <- diag(x)
  E[lower.tri(E, TRUE)] <- x[sequence(rev(seq_along(x)))]
  E
}

给了

> TIC(x)
     [,1] [,2] [,3] [,4]
[1,]    3    0    0    0
[2,]    4    3    0    0
[3,]    8    4    3    0
[4,]    9    8    4    3

和

【讨论】：

这很酷。我会注意到，在这个例子中，最大尺寸（100）的 slowest 方法的中位时间小于 0.01 秒 ...

【解决方案3】：

这是一个基本的 R 方式。

E <- diag(length(x))
apply(lower.tri(E, diag = TRUE), 2, function(i) {
  c(rep(0, nrow(E) - sum(i)), x)[seq_along(x)]
})
#     [,1] [,2] [,3] [,4]
#[1,]    3    0    0    0
#[2,]    4    3    0    0
#[3,]    8    4    3    0
#[4,]    9    8    4    3

性能测试

如果问题是关于更快的代码，这里是基准测试。

函数是我的和Ben Bolker's代码。

Rui <- function(x){
  E <- diag(length(x))
  inx <- seq_along(x)
  apply(lower.tri(E, diag = TRUE), 2, function(i) {
    c(rep(0, nrow(E) - sum(i)), x)[inx]
  })
}

Ben <- function(x){
  E <- matrix(0, nrow=length(x), ncol=length(x))
  diag(E) <- x[1]
  for (i in 2:length(x)) {
    E[row(E)==col(E)+i-1] <- x[i]
  }
  E
}

使用增加的向量大小进行测试并使用ggplot 绘图。

library(microbenchmark)
library(ggplot2)

test_speed <- function(n){
  out <- lapply(1:n, function(i){
    x <- sample(10*i)
    mb <- microbenchmark(
      Rui = Rui(x),
      Ben = Ben(x)
    )
    mb <- aggregate(time ~ expr, mb, median)
    mb$size <- 10*i
    mb
  })
  out <- do.call(rbind, out)
  out
}

res <- test_speed(10)

ggplot(res, aes(size, time, color = expr)) +
  geom_line() +
  geom_point() +
  scale_y_continuous(trans = "log10")

【讨论】：

最好将dapply 包含在您的函数中并对其进行测试。谢谢
时间单位是什么？我认为它们可能是纳秒，这可以将时间差异放在透视图上......（如果你要为真正的大矩阵执行此操作，你会想要使用稀疏矩阵，特别是三角形类之一...）
@BenBolker 是的，它们是纳秒级。

【解决方案4】：

这不是超级效率，但比您的解决方案更好。（效率低下的是我们每次都在构造row()/col() 矩阵并生成一个完整的逻辑矩阵，而不是做一些索引操作。）另一方面，对于length(x)==100（当我们达到 1000 时有点慢）。

E <- matrix(0, nrow=length(x), ncol=length(x))
diag(E) <- x[1]
for (i in 2:length(x)) {
   E[row(E)==col(E)+i-1] <- x[i]
}

可能有人编写了更高效的代码（在 Rcpp 中？）用于索引矩阵的子对角线/非对角线元素。

尽管速度很慢，但这个（IMO）的优点是它更容易理解；您还可以通过对行和列之间的关系提出不同的条件来将其调整为许多不同的模式。

【讨论】：