矩阵列表的高效子集化
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
循环,如果可能的话,像 mapply
和 lapply
这样的函数因为我必须调用这个函数数百万次。 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
我有一个这样的矩阵列表:
[[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 基础选择函数([[
、[
或 $
运算符)进行子集化?
如果您的矩阵都具有相同的维度,您可以将它们转换为 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