假设我有以下两个数组：
R <- 101
v <- array(0, dim <- c(R,2))
v[,1] <-runif(101)
t <- array(runif(5), dim <- c(5,2))
我想做的是为 v 第二列中的每个单元格分配以下函数的结果：
which.min(abs(v[r,1] - t[,1]))
因此，对于 v 的第二列中的每个单元格，我将有一个 1、2、3、4 或 5。我知道我可以在 v 的所有行 r 上使用 for 循环来做到这一点，但是有人知道一种向量化的方法这个操作让我不必求助于（相当慢的）for循环？
您可以扩展v和t：
V <- matrix(rep.int(v[,1],dim(t)[[1]]),ncol=dim(t)[[1]])
TT <- matrix(rep.int(t[,1],dim(v)[[1]]),ncol=dim(t)[[1]],byrow=T)
然后减去并取每一列的最大值：
max.col(-abs(V-TT))
尽管有名称，但并没有真正矢量化，如Vectorize调用lapply。但这给出了结果：
> Vectorize(function(r) which.min(abs(v[r,1] - t[,1])))(seq(nrow(v)))
## [1] 4 3 3 2 5 5 2 5 2 5 3 3 2 5 1 4 5 5 4 3 3 5 5 2 4 2 2 4 4 3 2 4 5 2
## [35] 2 3 2 4 4 1 5 5 2 3 2 4 5 5 3 5 2 4 4 2 4 5 5 5 5 5 4 3 3 5 5 3 2 3
## [69] 5 3 5 3 3 5 4 5 5 3 1 2 5 5 2 3 3 4 3 3 4 5 4 2 2 3 4 2 5 5 5 5 2
然后可以将此值分配给v[,2]。
我认为可以提供一种使用stepfun和结合的矢量化解决方案pmin，pmax所有这些都是矢量化的。它有点扭曲/复杂的逻辑，但值得付出所有努力。
使用stepfun++的优点pmin：pmax
首先，这个想法的灵感来自这里Jonathan Chang's的帖子。这里的小变化是您需要索引而不是差异。另外，我假设所有值都是正的（来自您的输入）。您可以将其扩展到具有负输入的向量，但如果需要，我将这项任务留给您。在我开始代码和基准测试之前，让我解释一下背后的想法是什么。runifstepfun
假设您有两个向量x（等价于v[,1]）和y（等价于t[,1]）。现在，让我们以这种方式排序并y创建一个：stepfunsorted y
y_sort <- sort(y)
step <- stepfun(y_sort, 0:length(y))
这对我们有多大帮助？查询step(a)会为您提供其中最大值的y_sort索引 < a。这可能需要一段时间才能理解。换句话说，值a位于step(a)和之间step(a) + 1的位置sorted y (y_sort)。现在，我们首先要弄清楚的是，这两个值中的哪一个最接近a. 这是通过提取索引step(a)和与这些索引对应step(a)+1的值并询问是否. 如果它是假的，那么，就是你的索引，反之亦然。其次，从from中取回原始索引，这可以通过使用in选项获取相应的排序索引来实现。y_sortabs(a-y_sort[step(a)]) > abs(a - y_sort[step(a)+1])step(a)yy_sortindex.return = TRUEsort
我同意以这种方式遵循这可能会非常复杂。但是检查代码并逐步运行它，并使用上面的文本跟随它（如有必要）。最好的部分是它a可以是一个向量，所以它非常快！现在开始代码。
# vectorised solution using stepfun
vectorise_fun1 <- function(x, y) {
y_sort <- sort(abs(y), index.return = TRUE)
y_sval <- y_sort$x
y_sidx <- y_sort$ix
# stepfun
step_fun <- stepfun(y_sval, 0:length(y))
ix1 <- pmax(1, step_fun(x))
ix2 <- pmin(length(y), 1 + ix1)
iy <- abs(x - y_sval[ix1]) > abs(x - y_sval[ix2])
# construct output
res <- rep(0, length(x))
res[iy] <- y_sidx[ix2[iy]]
res[!iy] <- y_sidx[ix1[!iy]]
res
}
# obtaining result
out_arun <- vectorise_fun1(v[,1], t[,1])
# (or) v[,2] <- vectorise_fun1(v[,1], t[,1])
# Are the results identical?
# Matthew's solution
vectorise_fun2 <- function(x, y) {
res <- Vectorize(function(r) which.min(abs(x[r] - y)))(seq(length(x)))
}
out_matthew <- vectorise_fun2(v[,1], t[,1])
# Jonathan's solution
vectorise_fun3 <- function(x, y) {
V <- matrix(rep.int(x, length(y)), ncol = length(y))
TT <- matrix(rep.int(y, length(x)), ncol = length(y), byrow = T)
max.col(-abs(V-TT))
}
out_jonathan <- vectorise_fun3(v[,1], t[,1])
# Are the results identical?
> all(out_arun == out_matthew)
[1] TRUE
> all(out_arun == out_jonathan)
[1] TRUE
那么，有什么意义呢？所有结果都是相同的，并且功能stepfun庞大且难以遵循。让我们取一个巨大的向量。
x <- runif(1e4)
y <- runif(1e3)
现在，让我们进行基准测试以查看优势：
require(rbenchmark)
> benchmark( out_arun <- vectorise_fun1(x,y),
out_matthew <- vectorise_fun2(x,y),
out_jonathan <- vectorise_fun3(x,y),
replications=1, order = "elapsed")
# test replications elapsed relative user.self
# 1 out_arun <- vectorise_fun1(x, y) 1 0.004 1.00 0.005
# 2 out_matthew <- vectorise_fun2(x, y) 1 0.221 55.25 0.169
# 3 out_jonathan <- vectorise_fun3(x, y) 1 1.381 345.25 0.873
# Are the results identical?
> all(out_arun == out_matthew)
[1] TRUE
> all(out_arun == out_jonathan)
[1] TRUE
因此，使用step_fun速度至少快 55 倍，最多快 345 倍！现在，让我们寻找更大的向量。
x <- runif(1e5)
y <- runif(1e4)
require(rbenchmark)
> benchmark( out_arun <- vectorise_fun1(x,y),
out_matthew <- vectorise_fun2(x,y),
replications=1, order = "elapsed")
# test replications elapsed relative user.self
# 1 out_arun <- vectorise_fun1(x, y) 1 0.052 1.000 0.043
# 2 out_matthew <- vectorise_fun2(x, y) 1 16.668 320.538 11.849
乔纳森的函数导致分配错误：
Error in rep.int(x, length(y)) :
cannot allocate vector of length 1000000000
而且这里的加速是320倍。