在 R 中计算 big.matrix 的行总和？

Posted 2023-04-18

【中文标题】在 R 中计算 big.matrix 的行总和？

【英文标题】：Computing row sums of a big.matrix in R?

【发布时间】：2014-07-10 22:32:24

【问题描述】：

我有一个大矩阵，大约有 6000 万行和 150 列（总共大约 90 亿个元素）。我已将此数据存储在big.matrix 对象中（来自包bigmemory）。现在，我希望计算每一行的总和，这是一个问题，因为big.matrix 是面向列的，所以据我所知，所有汇总函数都是面向列的（例如colsum、colmax 等.) 并且默认情况下没有可用于计算行总和的函数。我当然可以apply(x, 1, sum)，但这需要很长时间。我还可以逐列循环并使用矢量化添加来添加它们：

mysum <- rep(0, nrow(x))
for (i in seq(ncol(x))) 
  mysum <- mysum + x[,i]

但这仍然需要 20 多分钟，而且显然不是最理想的，因为它每次循环都会创建一个新的 6000 万元素向量。似乎必须有一些更快的方法来做到这一点。

编辑

我通过一次处理大约一百万行的块，并在这些块上调用 rowSums，然后连接结果，将这个时间缩短到了 10 分钟。不过，我仍然想知道是否有优化的方法来做到这一点。

【问题讨论】：

rowSums 对它不起作用吗？可以转置然后取colsum吗？

假设您有 numeric 数据，您指定的时间对应于大约 60 MB/s 的吞吐量（20 分钟内 72 GB 数据 = 3.6 GB 每分钟）。根据数据的存储位置，这可能非常接近物理极限。 读取该文件 (time cp file > /dev/null) 需要多长时间？

不是 R 的数字类型。这是一个big.matrix 的整数，所以我相信它在磁盘和内存中都存储得更紧凑。磁盘文件大约 30 GB，我不知道它是否在加载时将整个矩阵加载到内存中。您不能一次对整个事物进行操作，因为其中包含更多 .Machine$integer.max 元素。这就是为什么我将它放在big.matrix 中的原因。您不能快速转置big.matrix。就像我说的，数据结构是面向列的，所以转置它必须完全重建整个数据结构。

您可以随时修改Rcpp 库中的代码以执行rowSums 而不是colSums：gallery.rcpp.org/articles/using-bigmemory-with-rcpp

【参考方案1】：

为此，我编写了一些 C++ 代码，改编自 bigmemory Rcpp gallery:

rowSums.cpp

// [[Rcpp::depends(BH)]]
#include <Rcpp.h>
using namespace Rcpp;

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

#include <numeric>

// Logic for BigRowSums.
template <typename T>
NumericVector BigRowSums(XPtr<BigMatrix> pMat, MatrixAccessor<T> mat) 
    NumericVector rowSums(pMat->nrow(), 0.0);
    NumericVector value(1);
    for (int jj = 0; jj < pMat->ncol(); jj++) 
      for (int ii = 0; ii < pMat->nrow(); ii++) 
        value = mat[jj][ii];
        if (all(!is_na(value))) 
          rowSums[ii] += value[0];
           
         
       
    return rowSums;


// Dispatch function for BigRowSums
//
// [[Rcpp::export]]
NumericVector BigRowSums(SEXP pBigMat) 
    XPtr<BigMatrix> xpMat(pBigMat);

    switch(xpMat->matrix_type()) 
      case 1:
        return BigRowSums(xpMat, MatrixAccessor<char>(*xpMat));
      case 2:
        return BigRowSums(xpMat, MatrixAccessor<short>(*xpMat));
      case 4:
        return BigRowSums(xpMat, MatrixAccessor<int>(*xpMat));
      case 6:
        return BigRowSums(xpMat, MatrixAccessor<float>(*xpMat));
      case 8:
        return BigRowSums(xpMat, MatrixAccessor<double>(*xpMat));
      default:
        throw Rcpp::exception("unknown type detected for big.matrix object!");
       

在 R 中：

library(bigmemory)
library(Rcpp)
sourceCpp("rowSums.cpp")

m <- as.big.matrix(matrix(1:9, 3))
BigRowSums(m@address)
[1] 12 15 18

【讨论】：

这看起来很棒。如果它比我已经得到的更快（可能是），我会接受它作为答案。

我确实尝试创建一个文件备份的大矩阵来测试它，但它使我的计算机停止运行（它的创建），所以我杀死了它。我很想看看它的表现如何！

你不能使用 Eigen 或 Armadillo 吗？我不认为即使在 C++ 中做这样的嵌套循环也会非常快......