在R中索引矩阵的元素

Question

这个问题很傻，但我想知道我是否遗漏了什么。 假设有一个向量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]
}

Answer 1

这是一个使用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。

结果：

显然，性能很大程度上取决于矩阵的大小。对于大型矩阵，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

注意：我没有在比较中包含 George 的方法，因为从 Joseph 的结果来看，它似乎慢了很多。所以我的基准测试中比较的所有方法都只用基础 R 编写。

Answer 2

@docendodiscimus 回答的变体：您可以更改行和列索引，而不是转置，您可以通过将 lower.tri 包装在 which 中来获得：

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]

我猜如果矩阵很大，转置可能会很昂贵，这就是我考虑这个选项的原因。注意：约瑟夫伍德的回答表明我错了，因为转置方式在他的基准测试中更快。

---

（感谢@JosephWood：）您可以使用(n^2 - n)/2 + n，而不是使用sum(seq(n)) 进行枚举和求和。

Answer 3

library(miscTools)
k <- 1:15
triang(k, 5)

Answer 4

这是一个非常快速的基础 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 在 cmets 中指出的那样，t1 是一个整数矩阵，而其他的是数字矩阵。我应该知道这一点，因为 most important R functions 会从一开始就告诉我这些信息。

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 ...