从距离矩阵创建 dist 对象的内存高效方法答案

【问题标题】：Memory-efficient method to create dist object from distance matrix从距离矩阵创建 dist 对象的内存高效方法
【发布时间】：2020-05-11 09:44:46
【问题描述】：

我正在尝试从大距离矩阵创建dist 对象。我使用 stats::as.dist 时内存不足。例如，我在当前机器上有大约 128 Gb 可用，但 as.dist 在处理 73,000 x 73,000 矩阵（大约 42Gb）时内存不足。鉴于最终的dist 对象应该小于矩阵大小的一半（即它是下三角形，存储为向量），在我看来，应该可以以更节省内存的方式进行此计算方式 - 前提是我们注意不要创建大型中间对象，只需将输入的相关元素直接复制到输出即可。

查看getS3method('as.dist', 'default') 的代码，我发现它使用ans <- m[row(m) > col(m)] 进行计算，这需要创建与输入维数相同的row 和col 矩阵。

我认为我可以使用here 中的算法来改进这一点，以生成下三角形的索引。这是我使用这种方法的尝试。

as.dist.new = function(dm, diag = FALSE, upper = FALSE) {
  n = dim(dm)[1]
  stopifnot(is.matrix(dm))
  stopifnot(dim(dm)[2] == n)
  k = 1:((n^2 - n)/2)
  j <- floor(((2 * n + 1) - sqrt((2 * n - 1) ^ 2 - 8 * (k - 1))) / 2)
  i <- j + k - (2 * n - j) * (j - 1) / 2
  idx = cbind(i,j)
  remove(i,j,k)
  gc()
  d = dm[idx]

  class(d) <- "dist"
  attr(d, "Size") <- n
  attr(d, "call") <- match.call()
  attr(d, "Diag") <- diag
  attr(d, "Upper") <- upper
  d
}

这适用于较小的矩阵。这是一个简单的例子：

N = 10
dm = matrix(runif(N*N), N, N) 
diag(dm) = 0

x = as.dist(dm)
y = as.dist.new(dm)

但是，如果我们创建一个更大的距离矩阵，它会遇到与as.dist 相同的内存问题。

例如此版本在我的系统上崩溃：

N = 73000
dm = matrix(runif(N*N), N, N)
gc() 
diag(dm) = 0
gc()

as.dist.new(dm)

有没有人建议如何更有效地执行此操作？欢迎使用 R 或 Rcpp 解决方案。注意看this answer 到一个相关问题（从2 列位置数据生成全距离矩阵）似乎可以使用RcppArmadillo 来做到这一点，但我没有使用它的经验。

【问题讨论】：

请按照这里的建议尝试@mharinga 的解决方案：stackoverflow.com/questions/37407065/…
感谢@tushaR 的链接 - 该答案中的包可能有助于更快地生成距离矩阵。但我在做的是不同的。它将距离矩阵转换为 dist 对象。
这里没有惊喜。你有 42GB 矩阵 + 21GB 向量，已经是 63GB。您可以尝试使用例如将矩阵放在磁盘上来自包 {bigstatsr} 的 FBM（免责声明：我是作者）。
感谢@F.Privé 的建议 - 但我的问题不是磁盘空间不足（我当然可以增加交换大小），而是我应该已经有足够的内存可用矩阵加向量。因此，需要寻找更有效的解决方案来限制所需的磁盘交换量。
您没有足够的 RAM 用于 2 个对象 + 其他东西。这就是为什么我建议您将矩阵放在磁盘上。

标签： r rcpp

【解决方案1】：

嗯，遍历下三角的逻辑比较简单，如果你用 C++ 做，那么它也可以很快：

library(Rcpp)

sourceCpp(code='
// [[Rcpp::plugins(cpp11)]]

#include <cstddef> // size_t

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector as_dist(const NumericMatrix& mat) {
  std::size_t nrow = mat.nrow();
  std::size_t ncol = mat.ncol();
  std::size_t size = nrow * (nrow - 1) / 2;
  NumericVector ans(size);

  if (nrow > 1) {
    std::size_t k = 0;
    for (std::size_t j = 0; j < ncol; j++) {
      for (std::size_t i = j + 1; i < nrow; i++) {
        ans[k++] = mat(i,j);
      }
    }
  }

  ans.attr("class") = "dist";
  ans.attr("Size") = nrow;
  ans.attr("Diag") = false;
  ans.attr("Upper") = false;
  return ans;
}
')

as_dist(matrix(1:100, 10, 10))
    1  2  3  4  5  6  7  8  9
2   2                        
3   3 13                     
4   4 14 24                  
5   5 15 25 35               
6   6 16 26 36 46            
7   7 17 27 37 47 57         
8   8 18 28 38 48 58 68      
9   9 19 29 39 49 59 69 79   
10 10 20 30 40 50 60 70 80 90

【讨论】：

谢谢。适用于较小的矩阵。但是在较大的内存问题上遇到相同的内存问题。虽然我承认我不知道为什么 - 它看起来好像应该工作。我猜这与 R 创建传递或返回的对象的副本有关。但是我认为只有指针在 R 和 C++ 之间传递，所以我很困惑为什么会这样。
@dww 很奇怪，您收到的确切错误消息是什么？您使用的是哪个操作系统？
我使用系统监视器看到内存（RAM 和交换）逐渐填满，直到 R 会话在两者都达到 100% 时崩溃。在 64Gb Ram+ 64Gb Swap 上使用 Linux (Pop!OS)。
@dww 不知道哪里出了问题，也许可以尝试添加一些好的旧的Rcout << "foo" 调试语句，看看它是否真的运行代码或在分配ans 后立即中断？
我已经稍微改变了方法，现在正在使用一个函数（大致基于您的答案）直接从 lat 和 lon 向量计算 dist 对象，根本不生成中间距离矩阵。这似乎是最好的前进方式。稍后我会尝试您的调试建议，因为了解此版本失败的原因也很有趣。一种可能性——j 不应该从0 迭代到(ncol-1)？将循环运行到ncol 可能会超出范围。虽然我希望这会引发越界错误，而不仅仅是内存不足。

【解决方案2】：

我目前的解决方案是直接从 lat 和 lon 向量计算 dist 对象，根本不生成中间距离矩阵。在大型矩阵上，这比 geosphere::mdist() 后跟 stats::as.dist() 的“常规”管道快几百倍，并且只使用存储最终 dist 对象所需的内存。

以下 Rcpp 源代码基于使用 here 中的函数计算 c++ 中的半正弦距离，以及对 @Alexis method 的改编以遍历 c++ 中的下三角形元素。

#include <Rcpp.h>
using namespace Rcpp;

double distanceHaversine(double latf, double lonf, double latt, double lont, double tolerance){
  double d;
  double dlat = latt - latf;
  double dlon =  lont - lonf;
  d = (sin(dlat * 0.5) * sin(dlat * 0.5)) + (cos(latf) * cos(latt)) * (sin(dlon * 0.5) * sin(dlon * 0.5));
  if(d > 1 && d <= tolerance){
    d = 1;
  }
  return 2 * atan2(sqrt(d), sqrt(1 - d)) * 6378137.0;
}

double toRadians(double deg){
  return deg * 0.01745329251;  // PI / 180;
}

//-----------------------------------------------------------
// [[Rcpp::export]]
NumericVector calc_dist(Rcpp::NumericVector lat, 
                        Rcpp::NumericVector lon, 
                        double tolerance = 10000000000.0) {
  std::size_t nlat = lat.size();
  std::size_t nlon = lon.size();
  if (nlat != nlon) throw std::range_error("lat and lon different lengths");
  if (nlat < 2) throw std::range_error("Need at least 2 points");
  std::size_t size = nlat * (nlat - 1) / 2;
  NumericVector ans(size);
  std::size_t k = 0;
  double latf;
  double latt;
  double lonf;
  double lont;
  
  for (std::size_t j = 0; j < (nlat-1); j++) {
    for (std::size_t i = j + 1; i < nlat; i++) {
      latf = toRadians(lat[i]);
      lonf = toRadians(lon[i]);
      latt = toRadians(lat[j]);
      lont = toRadians(lon[j]);
      ans[k++] = distanceHaversine(latf, lonf, latt, lont, tolerance);
    }
  }
  
  return ans;
}

/*** R
as_dist = function(lat, lon, tolerance  = 10000000000.0) {
  dd = calc_dist(lat, lon, tolerance)
  attr(dd, "class") = "dist"
  attr(dd, "Size") = length(lat)
  attr(dd, "call") = match.call()
  attr(dd, "Diag") = FALSE
  attr(dd, "Upper") = FALSE
  return(dd)
}
*/

【讨论】：