查找大矩阵的行之间的最小距离：分配限制错误

Posted 2023-02-22

【中文标题】查找大矩阵的行之间的最小距离：分配限制错误

【英文标题】：Find lowest distances between rows of a large matrix: Allocation limit error

【发布时间】：2016-10-31 10:11:06

【问题描述】：

我想计算一个大矩阵的所有行之间的距离。对于每一行，我需要找到距离最短的另一行。最终输出应该是一个列表，其中包含距离最短的行的 ID（参见下面示例中的 low_dis_ids）。

我能够找到小样本的解决方案（下面的示例）。但是，我无法使用更大的样本量来执行这些步骤，因为具有距离的矩阵变得非常大。有没有办法省略这么大的矩阵？我只需要带有 ID 的列表（如 low_dis_ids）。

可重现的例子：

set.seed(123)

# Calculation of distances with small samplesize is working well
N <- 100
data_100 <- data.frame(x1 = rnorm(N, 5, 10),
                       x2 = rnorm(N, 5, 10),
                       x3 = rnorm(N, 5, 10),
                       x4 = rnorm(N, 5, 10),
                       x5 = rnorm(N, 5, 10))

# Matrix with all distances (no problem for the smaller samplesize)
dist_100 <- as.matrix(dist(data_100))

# Find the row with the smallest distance
for(i in 1:nrow(dist_100)) 
  dist_100[i, i] <- Inf


low_dis <- numeric()
for(i in 1:nrow(dist_100)) 
  low_dis[i] <- as.numeric(sort(dist_100[ , i]))[1]


low_dis_ids <- list()
for(i in 1:length(low_dis)) 
  low_dis_ids[[i]] <- as.numeric(names(dist_100[ , i][dist_100[ , i] == low_dis[i]]))


# low_dis_ids is the desired output and stores the rows with the smallest distances



# The same procedure is not working for larger samplesizes
N <- 100000
data_100000 <- data.frame(x1 = rnorm(N, 5, 10),
                          x2 = rnorm(N, 5, 10),
                          x3 = rnorm(N, 5, 10),
                          x4 = rnorm(N, 5, 10),
                          x5 = rnorm(N, 5, 10))
dist_100000 <- dist(data_100000)

# Error: cannot allocate vector of size 37.3 Gb

【参考方案1】：

您绝对可以避免由于使用dist 而创建的大型矩阵。一种这样的解决方案是一次计算一行的距离，我们可以编写一个函数，给定整个数据帧和行 id 找到哪一行对应于最小距离。例如：

f = function(rowid, whole)
  d = colSums((whole[rowid,] - t(whole))^2) # calculate distance
  d[rowid] = Inf # replace the zero
  which.min(d)

colSums 函数优化得相当好，因此相对较快。我怀疑which.min 也是一种更快且可能更简洁的方法来循环遍历距离向量。

为了制定一个适用于任何此类数据框的解决方案，我编写了另一个短函数，该函数将其应用于每一行并为您提供行 ID 向量

mindists = function(dat) do.call(c,lapply(1:nrow(dat),f,whole = as.matrix(dat)))

如果您想要列表而不是向量，只需省略 do.call 函数。我这样做是为了更容易检查输出是否符合您的预期。

all(do.call(c,low_dis_ids) == mindists(data_100))
[1] TRUE

这也适用于我笔记本电脑上的更大示例。它并不快，因为您正在对f 进行nrow(data) 调用，但它确实避免了创建一个大对象。那里可能有更好的解决方案，但这是第一个想到的有效解决方案。希望对您有所帮助。

编辑：

已编辑，因为 Roland 提供了一个额外的 C++ 答案 我对较小的数据集进行了快速基准测试。在这种情况下，C++ 的答案肯定更快。 这个答案的一些额外推销是我认为如果你是纯粹的 R 程序员（不需要学习 C++ 和 RCpp）更容易理解。 R 版本使用lapply 的并行替换来并行化是微不足道的。我会注意这并不是要从 Rolands 的回答中拿走，我个人喜欢 Rcpp，只是为了给这个问题的任何未来读者提供额外的信息。

【讨论】：

非常感谢杰米罗文！这正是我正在寻找的！就像你说的，Rolands 的代码比你的要快，但是因为我不知道如何使用 C++，所以我更喜欢你的解决方案。

@JoachimSchork，没问题。很高兴它有帮助。

【参考方案2】：

使用 Rcpp，因为基本 R 解决方案会太慢：

library(Rcpp)
library(inline)
cppFunction(
"  IntegerVector findLowestDist(const NumericMatrix X) 
    const int n = X.nrow();
    const int m = X.ncol();
    IntegerVector minind(n);
    NumericVector minsqdist(n);
    double d;
    for (int i = 0; i < n; ++i) 
      if (i == 0) 
        d = 0;
        for (int k = 0; k < m; ++k) 
          d += pow(X(i, k) - X(1, k), 2.0);

        
        minsqdist(i) = d;
        minind(i) = 1;
       else 
        d = 0;
        for (int k = 0; k < m; ++k) 
          d += pow(X(i, k) - X(0, k), 2.0);

        
        minsqdist(i) = d;
        minind(i) = 0;
      

      for (int j = 1; j < n; ++j) 
        if (i == j) continue;
        d = 0;
        for (int k = 0; k < m; ++k) 
          d += pow(X(i, k) - X(j, k), 2.0);

        
        if (d < minsqdist(i)) 
          minsqdist(i) = d;
          minind(i) = j;
        
      
    
    return minind + 1;
  "
)

all.equal(findLowestDist(as.matrix(data_100)),
          unlist(low_dis_ids))
#[1] TRUE

findLowestDist(as.matrix(data_100000))
#works

算法或许可以改进。

【讨论】：

非常感谢 Roland，我尝试了您的代码，它运行良好。令人惊讶的是这段代码的运行速度有多快。但是，由于我不知道如何使用 C++，我使用了 jamieRowens 解决方案，即使您的解决方案更快。

您可以通过跟踪最小距离来使 C++ 算法更快，并且仅在小于最小距离时继续添加到 d。只要它变大了，就不需要继续添加了。