gpt4 book ai didi

R:两个数值矩阵的快速笛卡尔积计算

转载 作者:行者123 更新时间:2023-12-04 12:00:12 29 4
gpt4 key购买 nike

我有两个大的数字矩阵,想在 R 中计算它们的笛卡尔积。有没有办法比我目前的方法具有更高的性能和更低的内存使用量?

编辑:我添加了一个 Rcpp 版本,它的性能已经比我第一个只有 R 的方法好得多。由于我对 Rcpp 或 RcppArmadillo 没有经验:是否有更快/更标准化的方法来编写此 Rcpp 函数?

m1 <- matrix(sample(0:9, size=110000, replace = TRUE), ncol = 110)
m2 <- matrix(sample(0:9, size=110000, replace = TRUE), ncol = 110)

#Current approach:
m3 <- apply(m1, 1, function(x) x * t(m2))
matrix(m3, ncol = 110, byrow = TRUE)

#EDIT - Rcpp approach
library(Rcpp)
#assuming ncol(m1) == ncol(m2)
cppFunction('IntegerMatrix cartProd(IntegerMatrix m1, IntegerMatrix m2) {
int nrow1 = m1.nrow(), ncol = m1.ncol(), nrow2 = m2.nrow();
int orow = 0;
IntegerMatrix out(nrow1 * nrow2, ncol);

for (int r1 = 0; r1 < nrow1; r1++) {
for (int r2 = 0; r2 < nrow2; r2++) {
for (int c = 0; c < ncol; c++){
out(orow, c) = m1(r1, c) * m2(r2, c);
}
orow++;
}
}
return out;
}')
m5 <- cartProd(m1, m2)

最佳答案

您所猜测的最佳方法是使用 C++ 来执行您想要的笛卡尔积。与纯 Rcpp 版本相比,尝试将代码移植到 Armadillo 会产生较小的速度提升,这明显快于编写的 R 版本。有关每种方法的效果的详细信息,请参阅最后的基准测试部分。

第一个版本几乎是直接移植到 Armadillo ,实际上比最初的纯 Rcpp 函数性能略差。第二个版本使用 Armadillo 内置 submatrix viewseach_row()利用就地评估的功能。为了实现与 Rcpp 版本的奇偶校验,请注意传递引用的使用和产生 const arma::imat& 的带符号整数类型的使用。 .这避免了两个大 integer 的深拷贝。匹配类型并建立引用的矩阵。

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

// --- Version 1

// [[Rcpp::export]]
arma::imat cartProd_arma(const arma::imat& m1, const arma::imat& m2) {
int nrow1 = m1.n_rows, ncol = m1.n_cols, nrow2 = m2.n_rows, orow = 0;
arma::imat out(nrow1 * nrow2, ncol);

for (int r1 = 0; r1 < nrow1; ++r1) {
for (int r2 = 0; r2 < nrow2; ++r2) {
out.row(orow) = m1.row(r1) % m2.row(r2);
orow++;
}
}
return out;
}

// --- Version 2

// [[Rcpp::export]]
arma::imat cartProd_arma2(const arma::imat& m1, const arma::imat& m2) {
int nrow1 = m1.n_rows, ncol = m1.n_cols, nrow2 = m2.n_rows, orow = 0;
arma::imat out(nrow1 * nrow2, ncol);

for (int r1 = 0; r1 < nrow1; ++r1) {
out.submat(orow, 0, orow + nrow2 - 1, ncol - 1) = m1.row(r1) % m2.each_row();
orow += nrow2;
}
return out;
}

快速验证与初始产品匹配的实现细节
all.equal( cartProd(m1, m2), cartProd_arma(m1, m2))
# [1] TRUE
all.equal( cartProd(m1, m2), cartProd_arma2(m1, m2))
# [1] TRUE

为了生成基准,我通过预先转置矩阵来稍微整理初始函数,以避免每次在每行调用 apply 时进行多次转置调用。此外,我还包含了@user20650 展示的功能技巧。
# OP's initial R only solution with slight modifications
op_R = function(m1, m2){
m2 <- t(m2)
m3 <- matrix(apply(m1, 1, function(x) x * m2), ncol = ncol(m1), byrow = TRUE)
}

# user20650's comment
so_comment <- function(m1, m2){
m4 <- matrix(rep(t(m1), each=nrow(m1)) * c(m2), ncol=nrow(m1))
}

因此,我们有以下微基准
library("microbenchmark")
out <- microbenchmark(op_r = op_R(m1, m2), so_comment_r = so_comment(m1, m2),
rcpp = cartProd(m1, m2), arma_v1 = cartProd_arma(m1, m2),
arma_v2 = cartProd_arma2(m1, m2),
times = 50)
out
# Unit: milliseconds
# expr min lq mean median uq max neval
# op_r 1615.6572 1693.0526 1793.0515 1771.7353 1886.0988 2053.7050 50
# so_comment_r 2778.0971 2856.6429 2950.5837 2936.7459 3021.4249 3344.4401 50
# rcpp 463.6743 482.3118 565.0525 582.1660 614.3714 699.3516 50
# arma_v1 597.9004 620.5888 713.4101 726.7572 783.4225 820.3770 50
# arma_v2 384.7205 401.9744 490.5118 503.5007 574.6840 622.9876 50

所以,由此我们可以看出 cartProd_arma2 ,子矩阵 Armadillo 实现,是最好的函数,紧随其后的是 cartProd ,纯 Rcpp 实现。

关于R:两个数值矩阵的快速笛卡尔积计算,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/41525236/

29 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com