矩阵列表的高效子集化

Efficient subsetting over a list of matrix

我有一个这样的矩阵列表:

[[1]]
     [,1] [,2] [,3]
[1,]  0.8  2.0  3.2
[2,]  2.0  3.2  4.4
[3,]  3.2  4.4  5.6

[[2]]
      [,1] [,2] [,3]
[1,] -1.95 1.00 3.95
[2,]  1.00 3.95 6.90
[3,]  3.95 6.90 9.85

[[3]]
      [,1] [,2] [,3]
[1,] -1.1  0.1  1.3
[2,]  0.1  1.3  2.5
[3,]  1.3  2.5  3.7

我想得到一个向量,它只包含列表中每个矩阵的特定行和列,行由向量选择。

我编写了一个代码来执行此操作:

xmin <- NULL
row <- c(2,3,2)
for(i in 1:3){
  xmin[i] <- lista[[i]][row[i] , 1]
}

问题是我想以最有效的方式进行选择,我想避免 for 循环,如果可能的话,像 mapplylapply 这样的函数因为我必须调用这个函数数百万次。 mapply 的计时略逊于 for 循环计时,但速度不够快。

有没有可能使用 R 基础选择函数([[[$ 运算符)进行子集化?

如果您的矩阵都具有相同的维度,您可以将它们转换为 3 维数组并使用索引。

lst <- list(matrix(1:9, 3, 3), matrix(10:18, 3, 3), matrix(19:27, 3, 3))
arr <- do.call(cbind, lst)
dim(arr) <- c(3, 3, 3)

getl <- function(row, col)
sapply(1:3, function(i) lst[[i]][row[i], col])

geta <- function(row, col)
arr[cbind(row, col, 1:3)]

> system.time(replicate(100000, getl(1, 2)))
   user  system elapsed 
   2.65    0.00    2.65 

> system.time(replicate(100000, geta(1, 2)))
    user  system elapsed 
   0.47    0.00    0.47 

我发现值得对一些替代品进行基准测试,包括 akrun 和 Hong Ooi。

1) 具有预分配结果的 "for" 循环:

ff1 = function(x, vec) 
{
    ans = numeric(length(x))

    for(i in seq_along(x)) ans[i] = x[[i]][vec[i], 1]

    return(ans)
}

2) 已编译的 "for" 循环:

cmpff1 = compiler::cmpfun(ff1)

3) A mapply:

ff2 = function(x, vec) mapply(function(elt, i) elt[i, 1], x, vec)

4) cbind使用矩阵索引 (Hong Ooi):

ffHO = function(x, vec) 
        "dim<-"(do.call(cbind, x), 
                c(dim(x[[1]]), length(x)))[cbind(vec, 1, seq_len(length(x)))] 

5) (4) 的修改,只添加一个属性:

ffHO2 = function(x, vec) 
       "dim<-"(unlist(x), 
               c(dim(x[[1]]), length(x)))[cbind(vec, 1, seq_len(length(x)))]

与比较:

myls = replicate(5e4, matrix(runif(100), 10, 10), simplify = FALSE)
vec = sample(1:10, 5e4, T)

ans1 = ff1(myls, vec)
ans2 = cmpff1(myls, vec)
ans3 = ff2(myls, vec)
ans4 = ffHO(myls, vec)
ans5 = ffHO2(myls, vec)

identical(ans1, ans2)
#[1] TRUE
identical(ans1, ans3)
#[1] TRUE
identical(ans1, ans4)
#[1] TRUE
identical(ans1, ans5)
#[1] TRUE

microbenchmark::microbenchmark(ff1(myls, vec), cmpff1(myls, vec), 
                               ff2(myls, vec), ffHO(myls, vec), 
                               ffHO2(myls, vec), times = 15)
#Unit: milliseconds
#              expr       min        lq    median        uq      max neval
#    ff1(myls, vec) 113.26685 132.36089 138.28047 147.97974 240.7101    15
# cmpff1(myls, vec)  51.23446  55.35398  58.18066  69.07220  82.4652    15
#    ff2(myls, vec) 119.44709 138.54739 145.66654 156.75227 219.9084    15
#   ffHO(myls, vec) 119.57063 130.52029 141.02867 149.21742 174.8242    15
#  ffHO2(myls, vec)  40.69163  41.31125  47.80939  48.55551 118.1069    15