时间戳

首页 > 解决方案 > 从距离矩阵创建 dist 对象的内存高效方法

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

问题描述

我正在尝试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输入相同维度的矩阵。

我想我也许可以使用这里的算法来改进这一点，以生成下三角形的索引。这是我使用这种方法的尝试。

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 解决方案。注意查看相关问题的答案（从 2 列位置数据生成完整距离矩阵）似乎可以使用 来执行此操作RcppArmadillo，但我没有使用它的经验。

标签： rrcpp

解决方案

好吧，遍历下三角的逻辑比较简单，如果你用 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

推荐阅读