使用big.matrix对象计算欧几里得距离矩阵

8
我在R中有一个类为big.matrix的对象,维度为778844 x 2。这些值都是整数(公里)。我的目标是使用big.matrix计算欧几里得距离矩阵,并将其作为big.matrix类的结果对象。我想知道是否有最佳方法来完成这个任务。
我选择使用big.matrix类的原因是内存限制。我可以将big.matrix转换为matrix类的对象,并使用dist()计算欧几里得距离矩阵。然而,dist()返回的对象大小无法分配到内存中。
编辑:
以下答案由bigmemory软件包的作者和维护者John W. Emerson提供:
您可以使用大型代数,但这也是使用sourceCpp()和非常简短易懂的Rcpp的好用例。但简而言之,我们甚至不尝试提供高级功能(除了我们实现的基本功能作为概念验证)。一旦开始讨论超出内存范围的大型问题,没有单个算法可以涵盖所有用例。
1个回答

9
这里有一种使用 RcppArmadillo 的方法。很多部分与 RcppGallery 示例 非常相似。这将返回一个带有相关成对(按行)欧几里得距离的big.matrix。我喜欢用一个包装器函数来包装我的big.matrix函数,以创建更清晰的语法(即避免@address和其他初始化)。
注意-由于我们使用的是bigmemory(因此关注RAM使用情况),我让这个示例返回了仅包含下三角元素的 N-1 x N-1 矩阵。你可以修改这个,但这就是我拼凑出来的东西。

euc_dist.cpp

// To enable the functionality provided by Armadillo's various macros,
// simply include them before you include the RcppArmadillo headers.
#define ARMA_NO_DEBUG

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo, BH, bigmemory)]]

using namespace Rcpp;
using namespace arma;

// The following header file provides the definitions for the BigMatrix
// object
#include <bigmemory/BigMatrix.h>

// C++11 plugin
// [[Rcpp::plugins(cpp11)]]

template <typename T>
void BigArmaEuclidean(const Mat<T>& inBigMat, Mat<T> outBigMat) {

  int W = inBigMat.n_rows;

  for(int i = 0; i < W - 1; i++){
      for(int j=i+1; j < W; j++){
          outBigMat(j-1,i) = sqrt(sum(pow((inBigMat.row(i) - inBigMat.row(j)),2)));
      }
  }
}

// [[Rcpp::export]]
void BigArmaEuc(SEXP pInBigMat, SEXP pOutBigMat) {
  // First we tell Rcpp that the object we've been given is an external
  // pointer.
  XPtr<BigMatrix> xpMat(pInBigMat);
  XPtr<BigMatrix> xpOutMat(pOutBigMat);


  int type = xpMat->matrix_type();
  switch(type) {
      case 1:
        BigArmaEuclidean(
            arma::Mat<char>((char *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
            arma::Mat<char>((char *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
        );
        return;

      case 2:
        BigArmaEuclidean(
          arma::Mat<short>((short *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
          arma::Mat<short>((short *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
        );
        return;

      case 4:
        BigArmaEuclidean(
          arma::Mat<int>((int *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
          arma::Mat<int>((int *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
        );
        return;

      case 8:
        BigArmaEuclidean(
          arma::Mat<double>((double *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
          arma::Mat<double>((double *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
        );
        return;

      default:
        // We should never get here, but it resolves compiler warnings.
        throw Rcpp::exception("Undefined type for provided big.matrix");
  }

}

我的小包装

bigMatrixEuc <- function(bigMat){
    zeros <- big.matrix(nrow = nrow(bigMat)-1,
                        ncol = nrow(bigMat)-1,
                        init = 0,
                        type = typeof(bigMat))
    BigArmaEuc(bigMat@address, zeros@address)
    return(zeros)
}

测试

library(Rcpp)
sourceCpp("euc_dist.cpp")

library(bigmemory)

set.seed(123)
mat <- matrix(rnorm(16), 4)
bm <- as.big.matrix(mat)

# Call new euclidean function
bm_out <- bigMatrixEuc(bm)[]

# pull out the matrix elements for out purposes
distMat <- as.matrix(dist(mat))
distMat[upper.tri(distMat, diag=TRUE)] <- 0
distMat <- distMat[2:4, 1:3]

# check if identical 
all.equal(bm_out, distMat, check.attributes = FALSE)
[1] TRUE

1
我运行了上述代码,并得到了一个名为 bm_out 的矩阵。在阅读包装函数时,我原本以为 bm_out 应该是一个 big.matrix。难道我错了吗?这个例子确实应该返回一个 matrix 吗?有没有办法让 bm_out 直接成为一个 big.matrix(而不是先得到一个 matrix 再通过 as.big.matrix 转换)? - Ricky
1
@Ricky 只需删除方括号 '[]',因为它们仅用于确认输出与“dist”相同。 - cdeterman

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接