rpy2 和标准 R 矩阵之间的差异

Question

我的目标

因此，我尝试使用 巨大 矩阵上的 rpy2 包从 python3 脚本中的 preprocessCore R 包 (R-3.2.1) 调用 normalize_quantiles 函数（ 10 GB 以上的文件）。我几乎有无限的记忆。当我通过 R 本身运行它时，我能够完成规范化并打印到输出：

require(preprocessCore);
all <- data.matrix(read.table("data_table.txt",sep="\t",header=TRUE));
all[,6:57]=normalize.quantiles(all[,6:57]);
write.table(all,"QN_data_table.txt",sep="\t",row.names=FALSE);

我正在尝试将它构建到一个 python 脚本中，该脚本还使用rpy2 python 包执行其他操作，但我在构建矩阵的方式上遇到了问题。下面是一个例子：

matrix = sample_list  # My 2-d python array containing the data.
v = robjects.FloatVector([ element for col in matrix for element in col ])
m = robjects.r['matrix'](v, ncol = len(matrix), byrow=False)
print("Performing quantile normalization.")
Rnormalized_matrix = preprocessCore.normalize_quantiles(m)
norm_matrix = np.array(Rnormalized_matrix)

return header, pos_list, norm_matrix

---

问题

这适用于较小的文件，但是当我在我的大文件上运行它时，它会因错误而死：rpy2.rinterface.RRuntimeError: Error: cannot allocate vector of size 9.7 Gb

我知道 R 的向量的最大大小是 8 Gb，这就解释了为什么会抛出上述错误。 rpy2 docs 说：

“矩阵是数组的一种特殊情况。与数组一样，必须记住这只是一个具有维度属性（行数、列数）的向量。”

我有点想知道它是如何严格遵守这一点的，所以我更改了我的代码以初始化我想要的大小的矩阵，然后遍历并将元素分配给值：

matrix = sample_list  # My 2-d python array of data.
m_count = 1
m = robjects.r['matrix'](0.0, ncol=len(matrix), nrow=len(matrix[0]))
for samp in matrix:
    i_count = 1
    for entry in samp:
        m.rx[i_count, m_count] = entry  # Assign the data to the element.
        i_count += 1
    m_count += 1

print("Performing quantile normalization.")

Rnormalized_matrix = preprocessCore.normalize_quantiles(m)
norm_matrix = np.array(Rnormalized_matrix)

return header, pos_list, norm_matrix

同样，这适用于较小的文件，但会因与之前相同的错误而崩溃。

所以我的问题是允许在 R 中分配巨大的矩阵但在 rpy 中导致问题的根本区别是什么？我需要用不同的方法来解决这个问题吗？我应该把它吸起来并在R中做吗？或者有没有办法绕过我遇到的问题？

Answer 1

从根本上说，R 是一种函数式语言。这意味着在R中做的时候

m[i, j] <- 123

正在发生的事情是这样的

assign_to_cell <- `[<-`
m <- assign_to_cell(m, i, j, 123)

参数是按值传递的。

这意味着应该返回一个新矩阵m，其中包含新值的(i,j) 单元格。为每个赋值复制m 的效率会相当低，尤其是当您遇到较大的对象时，因此 R 解释器有一个绝妙的技巧（有关详细信息，请参阅 R 的 C 源代码）：表达式的左侧与表达式的右侧进行比较，如果对象相同，解释器将知道复制是不必要的，并且可以“就地”完成修改。

现在有了rpy2，还有两点需要考虑：虽然 Python（大部分）通过引用传递参数，但嵌入式 R 引擎无法知道 Python 表达式左侧发生了什么。

表达式

m.rx[i_count, m_count] = entry 

忠实地构建了一个 R 调用，例如

m <- assign_to_cell(m, i, j, entry)

但是 R 在表达式左侧看不到前方的能力。结果是制作了m 的副本。每次修改。

但是，rpy2 可以轻松地将 R 中定义的向量、矩阵和数组移动到 Python 的传递引用世界。例如，这些 R 对象可以别名为对应的 numpy 对象（使用 asarra - 请参阅 http://rpy.sourceforge.net/rpy2/doc-2.0/html/numpy.html#low-level-interface)。记住 R 数组是按列优先顺序排列的，还可以计算索引并跳过别名到 numpy 数组并就地修改单元格：

m[idx] = entry 

注意：

我认为，如果我没记错的话，由于 R 索引是 32 位整数而导致的 8Gb 限制在几年前就被取消了。 你的无限记忆可能比你想象的要少。系统上的物理内存并不一定意味着可以将其全部分配在一个连续的块中。