矩阵列表上的有效子集答案

【问题标题】：Efficient subsetting over a list of matrix矩阵列表上的有效子集
【发布时间】：2015-02-27 12:59:32
【问题描述】：

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

[[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 循环，如果可能的话，像 mapply 和 lapply 这样的函数，因为我必须调用这个函数数百万次。 mapply 的时序略逊于for 循环时序，但速度不够快。

是否有任何可能的方法来使用 R 基础选择功能（[[、[ 或 $ 运算符）进行此子集化？

【问题讨论】：

您是否尝试过mapply(function(x,y) x[y,1] , lista, row) 或尝试行/列索引。 rC <- cbind(row,1); mapply([, lista,lapply(split(rC, row(rC)), t))
另一个选项是m1 <- do.call(rbind, lista); m1[cbind(row+seq(0,length(lista)*nrow(lista[[1]])-1, by=3),1)] #[1] 2.00 3.95 0.10

标签： r list matrix fuzzy-logic

【解决方案1】：

如果您的矩阵都具有相同的维度，您可以将它们转换为 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

【讨论】：

这个解决方案工作得很好，我所有的矩阵都有相同的维度。我有一个疑问，我在列表中的矩阵数量在算法执行期间总是相同的，但是这个数字是一个参数，我必须使用它，所以我必须为我的问题做：@987654322 @ 和 geta 函数中：geta <- function(row, col) arr[cbind(row, col, 1:n_matrices)]

【解决方案2】：

我发现值得对一些替代品进行基准测试，包括 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) 一个mapply：

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

4）cbinding 并使用矩阵索引（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

【讨论】：