如何利用稀疏矩阵优化对许多位移点的评估

问题描述

我目前正在尝试有效解决以下问题：

我有两个向量v1，v2以及向量化函数f。对于x中的每个v1，我想计算f(x - v2)的平均值。关于此问题的特殊之处在于f在许多输入上将返回零。

示例：

set.seed(0)

v1 <- rnorm(1000)
v2 <- rnorm(1000)

f <- function(x) {
  ret <- double(length(x)) + 1
  ret[abs(x) > 0.01] <- 0
  ret
}

solution_01 <- function(v1,v2,f) {
  ret <- numeric(length(v1))
  for (x in v2) {
    ret <- ret + f(v1 - x)
  }
  ret/length(v2)
}

solution_02 <- function(v1,f) {
  apply(matrix(f(outer(v1,`-`)),nrow=length(v1)),1,sum)/length(v2)
}

solution_03 <- function(v1,f) {
  rowSums(matrix(f(outer(v1,nrow=length(v1)))/length(v2)
}

solution_04 <- function(v1,f) {
  rowMeans(matrix(f(outer(v1,nrow=length(v1)))
}

s1 <- solution_01(v1,f)
s2 <- solution_02(v1,f)
s3 <- solution_03(v1,f)
s4 <- solution_04(v1,f)

all.equal(s1,s2)
all.equal(s2,s3)
all.equal(s3,s4)

bench::mark(solution_01(v1,f),solution_02(v1,solution_03(v1,solution_04(v1,f))

# Sparsity
eval_points <- outer(v1,`-`)
sum(f(eval_points) == 0)/length(eval_points)

如您所见，我已经实现了四种可能的解决方案。目前， naive 解决方案（使用for循环）是最快的。我认为这是因为其他实现依赖于outer，这需要花费一些时间来分配所需的内存。

如何优化此代码？有没有办法利用f(outer(v1,v_2))的稀疏性？

解决方法

使用排序向量可以有效地解决您提出的问题（我很欣赏这可能是对您实际尝试做的简化），尤其是当向量的长度超过1000个元素时。

由于您要在v2中的值的一小段距离内对v1中的值进行计数，因此您可以拥有一种算法，该算法前进到第一个向量，直到超出所考虑的元素范围在第二个中，然后切换到通过另一个向量前进。这样，您只需通过每个向量一次，而无需length(v1)次。

正如指出的那样，R在这种情况下效率不高，如果您希望它真的很快，则应在Rcpp中编写整个代码。

...

我想写一个有趣的算法，结果证明即使使用R编写，示例数据的速度也快将近10倍！

solution_05 <- function(v1,v2,f) {
  vs1 <- sort(v1)
  vs2 <- sort(v2)
  n1 <- integer(length(v1))
  i1 <- 1
  i2 <- 1
  l <- length(v1)
  for (x2 in vs2) {
    # advance i1 until v1 is in range below
    while (x2-vs1[i1] > 0.01 & i1 <= l) i1 <- i1+1
    if (i2 > i1) n1[i1:(i2-1)] <- n1[i1:(i2-1)]+1 else i2 <- i1
    # advance i2 until out of range adding 1 to n1[i2] each time
    while (vs1[i2]-x2 <= 0.01 & i2 <= l) {
      n1[i2] <- n1[i2]+1
      i2 <- i2+1
    }
  }
  s5 <- n1[rank(v1)]/length(v2)
  return(s5)
}