大数据方法来计算R中的距离集？

Posted 2023-02-14

【中文标题】大数据方法来计算R中的距离集？

【英文标题】：Big data ways to calculate sets of distances in R?

【发布时间】：2022-01-20 13:24:18

【问题描述】：

问题：我们需要一种大数据方法来计算点之间的距离。我们在下面概述了我们想用一个五观察数据框做的事情。但是，随着行数变大（> 100 万），这种特殊方法是不可行的。过去，我们使用 SAS 进行此类分析，但如果可能，我们更喜欢 R。 （注意：我不打算展示代码，因为虽然我在下面概述了一种在较小数据集上执行此操作的方法，但对于我们规模的数据来说，这基本上是一种不可能的方法。）

我们从商店的数据框开始，每个商店都有纬度和经度（虽然这不是空间文件，我们也不想使用空间文件）。

# you can think of x and y in this example as Cartesian coordinates
stores <- data.frame(id = 1:5,
                     x = c(1, 0, 1, 2, 0),
                     y = c(1, 2, 0, 2, 0))

stores
  id x y
1  1 1 1
2  2 0 2
3  3 1 0
4  4 2 2
5  5 0 0

对于每个商店，我们想知道 x 距离内的商店数量。在一个小的数据框中，这很简单。创建另一个所有坐标的数据框，重新合并，计算距离，如果距离小于 x，则创建一个指标，然后将指标相加（商店本身减去一个，距离为 0）。这将产生一个如下所示的数据集：

   id x y  s1.dist  s2.dist  s3.dist  s4.dist  s5.dist
1:  1 1 1 0.000000 1.414214 1.000000 1.414214 1.414214
2:  2 0 2 1.414214 0.000000 2.236068 2.000000 2.000000
3:  3 1 0 1.000000 2.236068 0.000000 2.236068 1.000000
4:  4 2 2 1.414214 2.000000 2.236068 0.000000 2.828427
5:  5 0 0 1.414214 2.000000 1.000000 2.828427 0.000000

当您（任意）将 1.45 以下视为“收盘价”时，您最终会得到如下所示的指标：

# don't include the store itself in the total
   id x y s1.close s2.close s3.close s4.close s5.close total.close
1:  1 1 1        1        1        1        1        1           4
2:  2 0 2        1        1        0        0        0           1
3:  3 1 0        1        0        1        0        1           2
4:  4 2 2        1        0        0        1        0           1
5:  5 0 0        1        0        1        0        1           2

最终产品应如下所示：

   id total.close
1:  1           4
2:  2           1
3:  3           2
4:  4           1
5:  5           2

感谢所有建议。

非常感谢

【问题讨论】：

数据集到底有多大？数据集是否太大而无法带入 R？有Hadoop 和其他分布式存储系统的数据存储选项。如果数据可以完全导入 R，有很多选择。您可以阅读其中一些选项here。

数据已经在 HPCC 上。问题是，要创建我上面描述的那种矩阵，它就像一个 1,000,000 x 1,000,000 数据帧，即使使用并行化和 HPC 也不理想。不过，如果我误解了你的建议，请纠正我。

我还应该补充一点，我们正在使用机密数据，因此我们可以使用或添加的软件包受到限制。任何连接到互联网的东西都是不允许的，如果我正确理解文档的话，这似乎排除了 Hadoop。

R 可以处理大约 2M 行（或列），因此您将不得不使用 HPC 中的聚类等方法。但是，由于您尝试执行的操作并不是特别复杂，您可能会发现data.table 包是您的最佳选择。我不确定您在坐标（即半正弦、文森蒂、欧几里得等）或比例（即英里、公里等）之间寻找什么度量标准，我只能提供一个包名字！

一百万点？太多了，您必须计算 n(n-1)/2 距离，即约 5000 亿距离

【参考方案1】：

你有什么理由不能循环而不是做一个大计算？

stores <- data.frame(id = 1:5,
                     x = c(1, 0, 1, 2, 0),
                     y = c(1, 2, 0, 2, 0))

# Here's a Euclidean distance metric, but you can drop anything you want in here
distfun <- function(x0, y0, x1, y1)
  sqrt((x1-x0)^2+(y1-y0)^2)


# Loop over each store
t(sapply(seq_len(nrow(stores)), function(i)
  distances <- distfun(x0 = stores$x[i], x1 = stores$x,
                       y0 = stores$y[i], y1 = stores$y)
  # Calculate number less than arbitrary cutoff, subtract one for self
  num_within <- sum(distances<1.45)-1
  c(stores$id[i], num_within)
))

生产：

     [,1] [,2]
[1,]    1    4
[2,]    2    1
[3,]    3    2
[4,]    4    1
[5,]    5    2

这将适用于您可以带入 R 的任何大小的数据集，但随着大小的增加它会变得更慢。这是在我的机器上运行几秒钟的 10,000 个条目的测试：

stores <- data.frame(id=1:10000, 
                     x=runif(10000, max = 10), 
                     y=runif(10000, max = 10))

          [,1] [,2]
    [1,]     1  679
    [2,]     2  698
    [3,]     3  618
    [4,]     4  434
    [5,]     5  402
...
 [9995,]  9995  529
 [9996,]  9996  626
 [9997,]  9997  649
 [9998,]  9998  514
 [9999,]  9999  667
[10000,] 10000  603

计算越多越慢（因为它必须在每对点之间运行，这总是 O(n^2)）但不知道您要计算的实际距离度量，我们无法优化再慢一点。

【讨论】：

这与以矢量化格式进行整个计算相同。您仍在重复计算。例如，一旦你计算了 1 和 2 之间的距离，你就再次计算 2 和 1 之间的距离，这使得这个函数的时间复杂度在 O(n^2) 中。而且我的朋友不会在超过 100 万行中工作

@Onyambu 是的，同意 - 但至少在 O(n^2) 的时间复杂度下它是可行的（可能一次创建数据库，而不是交互的东西？），而内存复杂度为O(n^2) 将需要 hardware 根本不存在 - 请参阅我对 jay 的回答的评论，估计 1M 行需要 ~4TB 的 RAM

此外，距离矩阵不能保证是对称的 - 在欧几里得空间中是对称的，但在许多研究领域中，A 和 B 之间的距离并不总是与B 和 A，如果是这种情况，就无法避免“重复”计算。

【参考方案2】：

您真的已经尝试过经典的dist() 功能了吗？核心是在C 中实现的，因此应该很快。

可能对矩阵的强制转换（无论如何都发生在 dist 中）已经花费了大量时间，也许可以立即将其作为矩阵而不是首先作为数据框读取。

M <- as.matrix(stores[-1])

dist(M, diag=TRUE, upper=TRUE)
#          1        2        3        4        5
# 1 0.000000 1.414214 1.000000 1.414214 1.414214
# 2 1.414214 0.000000 2.236068 2.000000 2.000000
# 3 1.000000 2.236068 0.000000 2.236068 1.000000
# 4 1.414214 2.000000 2.236068 0.000000 2.828427
# 5 1.414214 2.000000 1.000000 2.828427 0.000000

否则你可以试试这个C++ 实现，它基本上是@coatless's code 的副本。但是，我在 R 脚本中使用了 Rcpp 包。

library(Rcpp)
cppFunction('Rcpp::NumericMatrix calcPWD1 (const Rcpp::NumericMatrix & x)
  unsigned int outrows = x.nrow(), i = 0, j = 0;
  double d;
  Rcpp::NumericMatrix out(outrows,outrows);

  for (i = 0; i < outrows - 1; i++)
    Rcpp::NumericVector v1 = x.row(i);
    for (j = i + 1; j < outrows ; j ++)
      d = sqrt(sum(pow(v1-x.row(j), 2.0)));
      out(j,i)=d;
      out(i,j)=d;
    
  

  return out;
')

calcPWD1(M)
#          [,1]     [,2]     [,3]     [,4]     [,5]
# [1,] 0.000000 1.414214 1.000000 1.414214 1.414214
# [2,] 1.414214 0.000000 2.236068 2.000000 2.000000
# [3,] 1.000000 2.236068 0.000000 2.236068 1.000000
# [4,] 1.414214 2.000000 2.236068 0.000000 2.828427
# [5,] 1.414214 2.000000 1.000000 2.828427 0.000000

但是，基准测试显然支持dist，所以你应该试一试：

M_big <- M[sample(nrow(M), 1e4, replace=TRUE), ]  ## inflate to 10k rows
microbenchmark::microbenchmark(
  dist=dist(M_big, diag=TRUE, upper=TRUE),
  calcPWD1=calcPWD1(M_big),
  control=list(warmup=10L),
  times=3L
)
# Unit: milliseconds
#     expr       min        lq     mean   median        uq       max neval cld
#     dist  640.1861  660.1396  765.881  680.093  828.7284  977.3638     3  a 
# calcPWD1 1419.4106 1439.1353 1505.253 1458.860 1548.1736 1637.4873     3   b

请务必阅读@coatless's and Dirk Eddelbuettel's answers，他们在其中写了更多关于C、C++ 和R 的内容，并且还有其他版本的函数。

【讨论】：

dist 肯定会突破一百万个条目！在您的示例中使用 10k 条目运行它已经占用了大约 400MB 的内存，预计在 100k 时会增加到 40GB，并且在 OP 的 1M 行处需要 4TB 的内存。

我认为问题在于dist 是否会损坏或内存不足。