logo标签列表

在R中对矩阵元素进行索引

rmatrix
14

问题相当简单，但我想知道是否有什么遗漏。假设有一个包含一些数字的向量k，比如说：

> k
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

我想将此转换为矩阵

> m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    0    6    7    8    9
[3,]    0    0   10   11   12
[4,]    0    0    0   13   14
[5,]    0    0    0    0   15

我的第一个想法是使用upper.tri()，例如像这样m[upper.tri(m, diag = TRUE)] <- k，但那不会得到上面的矩阵。

有没有更加聪明的解决方案？下面是我的解决方案，但我不太自豪。

rows <- rep(1:5, 5:1)

cols1 <- rle(rows)$lengths


cols <- do.call(c, lapply(1:length(cols1), function(x) x:5))

for(i in 1:length(k)) {
  m[rows[i], cols[i]] <- k[i]
}
4个回答
15

这里有一种使用lower.tri和t转置结果的选项:

k <- 1:15
m <- matrix(0, 5,5)
m[lower.tri(m, diag = TRUE)] <- k
m <- t(m)
m 
#     [,1] [,2] [,3] [,4] [,5]
#[1,]    1    2    3    4    5
#[2,]    0    6    7    8    9
#[3,]    0    0   10   11   12
#[4,]    0    0    0   13   14
#[5,]    0    0    0    0   15

微基准测试

由于Joseph的基准测试存在一些混淆，这里提供另一个测试。我对大小为10*10、100*100、1000*1000、10000*10000的矩阵进行了测试。

结果:

pic

显然，性能严重依赖于矩阵的大小。对于大矩阵，Joseph的方法表现最快，而对于小矩阵，我的方法是最快的。请注意，这并未考虑内存效率。

可复制的基准测试:

Joseph <- function(k, n) {
  y <- 1L
  t <- rep(0L,n)
  j <- c(y, sapply(1:(n-1L), function(x) y <<- y+(n+1L)-x))
  t(vapply(1:n, function(x) c(rep(0L,x-1L),k[j[x]:(j[x]+n-x)]), t, USE.NAMES = FALSE))
}

Frank <- function(k, n) {
  m = matrix(0L, n, n)
  m[ which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1] ] = k
  m
}

docendo <- function(k,n) {
  m <- matrix(0L, n, n)
  m[lower.tri(m, diag = TRUE)] <- k
  t(m)
}

library(microbenchmark)
library(data.table)
library(ggplot2)
n <- c(10L, 100L, 1000L, 10000L)
k <- lapply(n, function(x) seq.int((x^2 + x)/2))

b <- lapply(seq_along(n), function(i) {
  bm <- microbenchmark(Joseph(k[[i]], n[i]), Frank(k[[i]], n[i]), docendo(k[[i]], n[i]), times = 10L)
  bm$n <- n[i]
  bm
})

b1 <- rbindlist(b)

ggplot(b1, aes(expr, time)) +
  geom_violin() +
  facet_wrap(~ n, scales = "free_y") +
  ggtitle("Benchmark for n = c(10L, 100L, 1000L, 10000L)")

检查结果的相等性：

all.equal(Joseph(k[[1]], n[1]), Frank(k[[1]], n[1]))
#[1] TRUE
all.equal(Joseph(k[[1]], n[1]), docendo(k[[1]], n[1]))
#[1] TRUE

注意：由于Joseph的结果显示，乔治的方法似乎要慢得多，因此我没有将其方法包括在比较中。因此，在我的基准测试中，所有比较的方法都仅使用基本R编写。

我知道一定有办法可以使用upper.tri或lower.tri来实现。谢谢！ - Theodor
感谢您提供详尽的基准测试。图表也是一个非常好的补充！您和@Frank提供的解决方案更加直观（在我看来更加简洁），对于大多数情况可能更有用。 - Joseph Wood
1
@JosephWood，这是一个有趣的问题，我对你的解决方案在处理大矩阵时表现得如此出色感到有些惊讶。希望你能获得更多的赞同。 - talat
@docendo，我不太担心投票。我很享受想出解决方案的过程，而且是的，我也很惊讶它运行得如此之快。在第一次对一个非常大的示例进行基准测试后，我确信我意外地测试了一个较小的矩阵，因为我使用了不同的k索引（即k1，k2等）。基本R总是让我感到惊奇。 - Joseph Wood
11
< p >一种@docendodiscimus答案的变种：不必转置，可以改变行和列的索引，这可以通过将lower.tri包装在which中获得: < /p >
n = 5
m = matrix(0, n, n)

m[ which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1] ] = seq(sum(seq(n)))


     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    0    6    7    8    9
[3,]    0    0   10   11   12
[4,]    0    0    0   13   14
[5,]    0    0    0    0   15

要了解它的工作原理，请看左侧步骤：

  • lower.tri(m, diag=TRUE)
  • which(lower.tri(m, diag=TRUE), arr.ind=TRUE)
  • which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1]

如果矩阵很大，转置可能会很耗时，这就是为什么我会考虑这个选项。注意：Joseph Wood的答案表明我是错误的，因为在他的基准测试中，转置的方式更快。

(感谢 @JosephWood:) 不要用 sum(seq(n)) 枚举并求和，你可以使用 (n^2 - n)/2 + n。

1
注意：(n^2 - n)/2 + n 等于 sum(seq(n))。很好的答案！ - Joseph Wood
8
library(miscTools)
k <- 1:15
triang(k, 5)
6

这里是一个非常快速的基础R解决方案：

更新

我稍微修改了代码，只调用了一次vapply而不是之前的sapply/vapply组合（我也取消了USE.NAMES=FALSE，因为它似乎没有任何区别）。虽然这样更加简洁，但在我的机器上并没有显著改变时间（我用图表重新运行了docendo的基准测试，结果几乎相同）。

Triangle1 <- function(k,n) {
    y <- -n
    r <- rep(0L,n)
    t(vapply(1:n, function(x) {y <<- y+n+2L-x; c(rep(0L,x-1L),k[y:(y+n-x)])}, r))
}

以下是一些时间安排：

Triangle2 <- function(k,n) {
    m <- matrix(0, n,n)
    m[lower.tri(m, diag = TRUE)] <- k
    t(m)
}

Triangle3 <- function(k, n) {
    m = matrix(0, n, n)
    m[ which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1] ] = k   ## seq(sum(seq(n)))  for benchmarking
    m
}

k2 <- 1:50005000
n2 <- 10^4

system.time(t1 <- Triangle1(k2,n2))
user  system elapsed           ## previously   user  system elapsed
2.29    0.08    2.41           ##              2.37    0.13    2.52

system.time(t2 <- Triangle2(k2,n2))
user  system elapsed 
5.40    0.91    6.30

system.time(t3 <- Triangle3(k2,n2))
user  system elapsed 
7.70    1.03    8.77 

system.time(t4 <- triang(k2,n2))
user  system elapsed 
433.45    0.20  434.88

有一件事让我感到困惑的是，Triangle1 产生的对象大小只有其他解决方案的一半。

object.size(t1)
400000200 bytes

object.size(t2)   ## it's the same for t3 and t4
800000200 bytes

当我进行一些检查时，情况只会更加混乱。

all(sapply(1:ncol(t1), function(x) all(t1[,x]==t2[,x])))
[1] TRUE

class(t1)
[1] "matrix"
class(t2)
[1] "matrix"

attributes(t1)
$dim
[1] 10000 10000
attributes(t2)
$dim
[1] 10000 10000

## not sure what's going on here
identical(t1,t2)
[1] FALSE

identical(t2,t3)
[1] TRUE

如评论中@Frank所指出的，t1是一个整数矩阵，而其他的是数值型。我应该知道这一点，因为其中一个最重要的R函数会从一开始就告诉我这个信息。
str(t1)
int [1:10000, 1:10000] 1 0 0 0 0 0 0 0 0 0 ...
str(t2)
num [1:10000, 1:10000] 1 0 0 0 0 0 0 0 0 0 ...
1
“Triangle3” 计算 k，而其他的则是免费给出的，对吧？我不知道这需要多长时间，但它使基准测试看起来不对称。 - Frank
1
@Frank，抱歉。现在的时间反映了k被传递到Triangle3。我要说seq(sum(seq(n)))非常快。对于n = 10^4，它在我的机器上注册为0.06。 - Joseph Wood
2
t1比较小且不相同，因为它是一个整数矩阵。 - Frank
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接