除了使用 R 基函数之外，是否有一种有效的方法来获得 "pmax"？
Is there an efficient way to obtain "pmax" other than using the R base function?
我想使用 Rcpp 创建一个函数，它可以胜过 R base 中的 pmax 函数。
我还尝试处理 Rcpp 函数内的缺失值，这可能不是一个好主意。
所有向量都必须有一些缺失值，并且它们都是正值。这就是我将 missing 重新编码为 -1 的原因，因此我可以将其添加回去，以防万一所有值都缺失时最大值不存在。
这是我的第一次尝试，但还没有成功：
library("benchr")
library("Rcpp")
Pmax <- function(...) {
argd_list <- list(...)
cppFunction("
NumericVector cpp_pmax(List args) {
List args0 = args[0];
int n_arg = args.length();
int n_vec = args0.length();
NumericVector out(n_vec);
out = args[0];
for (int i = 1; i < n_arg; ++i) {
NumericVector pa(n_vec);
pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (R_IsNA(out[j])) {
out[j] = -1;
}
if (R_IsNA(pa[j])) {
pa[j] = -1;
}
out[j] = std::max(out[j], pa[j]);
}
}
for (int j = 0; j < n_vec; ++j) {
if (out[j] == -1) {
out[j] = NA_REAL;
}
}
return out;
}
")
output <- cpp_pmax(argd_list)
return(output)
}
n <- 200000
x1 <- sample(0:1, n, replace = TRUE)
y1 <- sample(0:1, n, replace = TRUE)
z1 <- sample(0:1, n, replace = TRUE)
x1[sample(1:n, 90)]<-NA
y1[sample(1:n, 60)]<-NA
z1[sample(1:n, 70)]<-NA
pm1 <- Pmax(x1, y1, z1)
pm2 <- pmax(x1, y1, z1, na.rm = TRUE)
all(pm1 == pm2)
benchr::benchmark(pmax(x1, y1, z1, na.rm = TRUE),
Pmax(x1, y1, z1))
Benchmark summary:
Time units : milliseconds
expr n.eval min lw.qu median mean up.qu max total relative
pmax(x1, y1, z1, na.rm = TRUE) 100 1.34 1.37 1.39 1.44 1.46 1.74 144 1.00
Pmax(x1, y1, z1) 100 13.30 13.50 13.80 19.90 15.70 67.50 1990 9.88
编辑：
我删除了一些循环，只是在 Rcpp 之外用 NA 替换了 -1，它加快了一点，但仍然没有超过 R base pmax。
虽然 Rcpp::pmax 是一个很好的实现，但它只处理两个向量，不确定它是否可以处理缺失值。当缺少值时，我得到了不同的结果。
第二次尝试是：
Pmax1 <- function(...) {
args_list <- list(...)
cppFunction("
NumericVector cpp_pmax(List args) {
List args0 = args[0];
int n_arg = args.length();
int n_vec = args0.length();
NumericVector out(n_vec);
out = args[0];
for (int i = 1; i < n_arg; ++i) {
NumericVector pa(n_vec);
pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (R_IsNA(out[j])) {
out[j] = -1;
}
if (R_IsNA(pa[j])) {
pa[j] = -1;
}
out[j] = std::max(out[j], pa[j]);
}
}
return out;
}
")
output <- cpp_pmax(args_list)
output[output == -1] <- NA
return(output)
}
Pmax2 <- function(...) {
args_list <- list(...)
cppFunction("
NumericVector cpp_pmax(List args) {
NumericVector out = args[0];
int n_arg = args.length();
int n_vec = out.length();
for (int j = 0; j < n_vec; ++j) {
if (NumericVector::is_na(out[j])) out[j] = -1;
}
for (int i = 1; i < n_arg; ++i) {
NumericVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (NumericVector::is_na(pa[j])) pa[j] = -1;
out[j] = std::max(out[j], pa[j]);
}
}
return out;
}
")
output <- cpp_pmax(args_list)
output[output == -1] <- NA
return(output)
}
n <- 200000
x <- sample(0:5, n, replace = TRUE)
y <- sample(0:5, n, replace = TRUE)
z <- sample(0:5, n, replace = TRUE)
w <- sample(0:5, n, replace = TRUE)
x[sample(1:n, 900)]<-NA
y[sample(1:n, 600)]<-NA
z[sample(1:n, 700)]<-NA
z[sample(1:n, 800)]<-NA
benchr::benchmark(pmax(x, y, z, w, na.rm = TRUE),
Pmax1(x, y, z, w),
Pmax2(x, y, z, w))
Benchmark summary:
Time units : milliseconds
expr n.eval min lw.qu median mean up.qu max total relative
pmax(x, y, z, w, na.rm = TRUE) 100 2.38 2.43 2.46 2.46 2.48 2.6 246 1.00
Pmax1(x, y, z, w) 100 16.00 16.90 17.20 19.40 17.70 61.2 1940 6.98
Pmax2(x, y, z, w) 100 9.44 9.74 9.90 11.30 10.10 45.6 1130 4.02
有没有人知道如何让它比 R base pmax 更快？
我们的想法是拥有一个通用函数来处理不同数量的向量，所有这些都在 Rcpp 函数中。
根据@DirkEddelbuettel 和@Cole 的回答更新
感谢您帮助优化代码。受@DirkEddelbuettel 和@Cole 回答的启发，我只是添加 Rcpp::pmax 来删除其中一个循环，它也有助于加快速度。
library("bench")
library("Rcpp")
cppFunction("
IntegerVector cpp_pmax1(List args) {
IntegerVector tmp = args[0];
IntegerVector out = clone(tmp);
int n_arg = args.length();
int n_vec = out.length();
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (pa[j] > out[j]) out[j] = pa[j];
}
}
return out;
}
")
cppFunction("
IntegerVector cpp_pmax2(List args) {
IntegerVector tmp = args[0];
IntegerVector out = clone(tmp);
int n_arg = args.length();
int n_vec = out.length();
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
out = pmax(out, pa);
}
return out;
}
")
Pmax1 <- function(...) {
cpp_pmax1(list(...))
}
Pmax2 <- function(...) {
cpp_pmax2(list(...))
}
n <- 200000
x <- sample(0:5, n, replace = TRUE)
y <- sample(0:5, n, replace = TRUE)
z <- sample(0:5, n, replace = TRUE)
w <- sample(0:5, n, replace = TRUE)
k <- sample(0:5, n, replace = TRUE)
x[sample(1:n, 900)] <- NA
y[sample(1:n, 600)] <- NA
z[sample(1:n, 700)] <- NA
w[sample(1:n, 800)] <- NA
k[sample(1:n, 800)] <- NA
pm0 <- pmax(x, y, z, w, k, na.rm = TRUE)
pm1 <- Pmax1(x, y, z, w, k)
pm2 <- Pmax2(x, y, z, w, k)
benchr::benchmark(pmax(x, y, z, w, k, na.rm = TRUE),
Pmax1(x, y, z, w, k),
Pmax2(x, y, z, w, k))
Benchmark summary:
Time units : microseconds
expr n.eval min lw.qu median mean up.qu max total relative
pmax(x, y, z, w, k, na.rm = TRUE) 100 2880 2900 2920 3050 3080 8870 305000 5.10
Pmax1(x, y, z, w, k) 100 2150 2180 2200 2310 2350 8060 231000 3.85
Pmax2(x, y, z, w, k) 100 527 558 572 812 719 7870 81200 1.00
谢谢！
顺便说一句，请注意 Rcpp 糖已经有 Rcpp::pmax():
> library(Rcpp)
> cppFunction("NumericVector pm(NumericVector x, NumericVector y) {
+ return pmax(x,y);}")
> pm(10.0*(1:10), rep(50, 10))
[1] 50 50 50 50 50 60 70 80 90 100
> pm(10.0*(1:10), c(rep(50, 8), NA, 50))
[1] 50 50 50 50 50 60 70 80 NA 100
>
可能还有另一个更通用的功能的范围，但希望这也可以作为基准对您有所帮助。
编辑： 在我的第一个版本中，当我打算调用 pm()（使用 Rcpp::pmax()）时，我不小心调用了 pmax()。结果一样。
pm() 和 pmax() 与人们预期的速度顺序大致相同，因为两者都是矢量化的：
> library(microbenchmark)
> set.seed(123)
> x <- cumsum(rnorm(1e6))
> y <- cumsum(rnorm(1e6))
> microbenchmark(pmax(x,y), pm(x,y))
Unit: milliseconds
expr min lq mean median uq max neval cld
pmax(x, y) 3.94342 4.07488 4.66378 4.15433 5.39961 7.81931 100 a
pm(x, y) 3.58781 3.68886 4.74249 3.75815 5.38444 22.31268 100 a
>
我想你可以尝试 fcoalesce + fifelse（都来自 data.table 包）来定义你的 Pmax 函数，如下所示
Pmax <- function(..., na.rm = FALSE) {
u <- list(...)
if (na.rm) {
return(
Reduce(function(x, y) {
x <- fcoalesce(x, y)
y <- fcoalesce(y, x)
fifelse(x <= y, y, x)
}, u)
)
}
Reduce(function(x, y) fifelse(x <= y, y, x), u)
}
Benchmark（使用 OP post 中的数据进行测试）
- 如果启用
na.rm = TRUE，Pmax 比基本 R pmax 稍慢
> microbenchmark::microbenchmark(
+ pmax(x1, y1, z1, na.rm = TRUE),
+ Pmax(x1, y1, z1, na.rm = TRUE),
+ check = "equivalent",
+ unit = "relati ..." ... [TRUNCATED]
Unit: relative
expr min lq mean median uq
pmax(x1, y1, z1, na.rm = TRUE) 1.000000 1.00000 1.000000 1.000000 1.000000
Pmax(x1, y1, z1, na.rm = TRUE) 1.428545 1.87539 1.974959 2.022579 2.094833
max neval
1.000000 100
1.387139 100
- 如果您使用默认的
na.rm 选项，您会发现 Pmax 比基础 R pmax 稍快
> microbenchmark::microbenchmark(
+ pmax(x1, y1, z1),
+ Pmax(x1, y1, z1),
+ check = "equivalent",
+ unit = "relative"
+ )
Unit: relative
expr min lq mean median uq max neval
pmax(x1, y1, z1) 1.387953 1.32482 1.053983 1.220124 1.143867 0.266205 100
Pmax(x1, y1, z1) 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 100
从 bench::mark 中可以看出内存分配似乎存在一些问题。
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2(x, y, z, w))
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc
## <bch:expr> <bch:t> <bch:t> <dbl> <bch:byt>
## 1 pmax(x, y, z, w, na.rm = TRUE) 5.79ms 6.28ms 157. 781.3KB
## 2 Pmax2(x, y, z, w) 39.56ms 54.48ms 19.7 9.18MB
内存强制
与基数 pmax() 相比，内存分配 是 10 倍。您的 rcpp 比较直接，因此这暗示存在某种强制。在查看示例数据时，您正在将整数向量发送到数字签名。这造成了代价高昂的胁迫。让我们更新签名和代码以期待 IntegerVectors。为此，我只是将所有内容从 NumericVector 更改为 IntegerVector。
expression min median `itr/sec` mem_alloc
<bch:expr> <bch:t> <bch:t> <dbl> <bch:byt>
1 pmax(x, y, z, w, na.rm = TRUE) 1.89ms 2.33ms 438. 781.3KB
2 Pmax2_int(x, y, z, w) 37.42ms 49.88ms 17.6 2.32MB
重新编译
OP 代码在较大的函数代码中包含 cppFunction。除非我们需要在每个循环中重新编译它，否则我们可以编译然后从 R 中调用编译后的代码。这是此数据集大小的最大性能提升。
cppFunction("
IntegerVector cpp_pmax_pre(List args) {
IntegerVector out = args[0];
int n_arg = args.length();
int n_vec = out.length();
for (int j = 0; j < n_vec; ++j) {
if (IntegerVector::is_na(out[j])) out[j] = -1;
}
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (IntegerVector::is_na(pa[j])) pa[j] = -1;
out[j] = std::max(out[j], pa[j]);
}
}
return out;
}
")
Pmax2_int_pre <- function(...) {
args_list <- list(...)
output <- cpp_pmax_pre(args_list)
output[output == -1] <- NA
return(output)
}
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2_int_pre(x, y, z, w))
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc
## <bch:expr> <bch:> <bch:> <dbl> <bch:byt>
## 1 pmax(x, y, z, w, na.rm = TRUE) 2.31ms 2.42ms 397. 781.3KB
## 2 Pmax2_int_pre(x, y, z, w) 2.48ms 3.55ms 270. 2.29MB
更多内存和小优化
最后，我们还有更多的内存分配。这暗示我们可以做更多 - 在这种情况下，我们应该在 rcpp 中更新 NA_REAL。相关的，我们可以优化一下循环赋值。
cppFunction("
IntegerVector cpp_pmax_final(List args) {
IntegerVector out = args[0];
int n_arg = args.length();
int n_vec = out.length();
for (int j = 0; j < n_vec; ++j) {
if (IntegerVector::is_na(out[j])) out[j] = -1;
}
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
// simplify logic; if the element is not na and is greater than the out, update out.
if (!IntegerVector::is_na(pa[j]) & pa[j] > out[j]) out[j] = pa[j];
}
}
// update now in Rcpp instead of allocating vectors in R
for (int i = 0; i < n_vec; i++) {
if(out[i] == -1) out[i] = NA_INTEGER;
}
return out;
}
")
Pmax2_final <- function(...) {
cpp_pmax_final(list(...))
}
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2_final(x, y, z, w))
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc
## <bch:expr> <bch:> <bch:> <dbl> <bch:byt>
## 1 pmax(x, y, z, w, na.rm = TRUE) 2ms 2.08ms 460. 781.3KB
## 2 Pmax2_final(x, y, z, w) 1.19ms 1.45ms 671. 2.49KB
我们做到了*！我确信可能会有一些小的优化 - 我们访问了 pa[j] 三次，因此可能值得分配给一个变量。
奖金 - NA_INTEGER
根据Rcpp for Everyone，NA_INTEGER 应该等于最小整数值-2147483648。使用这个，我们可以删除 NA 的替换，因为 我们可以在处理 int 数据类型时直接与 NA 进行比较。
在实现过程中，我还发现了前一部分的一个问题——我们需要克隆初始参数，这样我们就不会不小心通过引用更改它。不过，我们仍然比基础 pmax().
稍微快一点
cppFunction("
IntegerVector cpp_pmax_last(List args) {
IntegerVector tmp = args[0];
IntegerVector out = clone(tmp);
int n_arg = args.length();
int n_vec = out.length();
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (pa[j] > out[j]) out[j] = pa[j];
}
}
return out;
}
")
Pmax2_last <- function(...) {
cpp_pmax_last(list(...))
}
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2_last(x, y, z, w),
)
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc `gc/sec`
## <bch:expr> <bch:> <bch:> <dbl> <bch:byt> <dbl>
## 1 pmax(x, y, z, w, na.rm = TRUE) 5.98ms 6.36ms 154. 781KB 0
## 2 Pmax2_last(x, y, z, w) 5.09ms 5.46ms 177. 784KB 0
我想使用 Rcpp 创建一个函数，它可以胜过 R base 中的 pmax 函数。 我还尝试处理 Rcpp 函数内的缺失值，这可能不是一个好主意。 所有向量都必须有一些缺失值，并且它们都是正值。这就是我将 missing 重新编码为 -1 的原因，因此我可以将其添加回去，以防万一所有值都缺失时最大值不存在。
这是我的第一次尝试，但还没有成功：
library("benchr")
library("Rcpp")
Pmax <- function(...) {
argd_list <- list(...)
cppFunction("
NumericVector cpp_pmax(List args) {
List args0 = args[0];
int n_arg = args.length();
int n_vec = args0.length();
NumericVector out(n_vec);
out = args[0];
for (int i = 1; i < n_arg; ++i) {
NumericVector pa(n_vec);
pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (R_IsNA(out[j])) {
out[j] = -1;
}
if (R_IsNA(pa[j])) {
pa[j] = -1;
}
out[j] = std::max(out[j], pa[j]);
}
}
for (int j = 0; j < n_vec; ++j) {
if (out[j] == -1) {
out[j] = NA_REAL;
}
}
return out;
}
")
output <- cpp_pmax(argd_list)
return(output)
}
n <- 200000
x1 <- sample(0:1, n, replace = TRUE)
y1 <- sample(0:1, n, replace = TRUE)
z1 <- sample(0:1, n, replace = TRUE)
x1[sample(1:n, 90)]<-NA
y1[sample(1:n, 60)]<-NA
z1[sample(1:n, 70)]<-NA
pm1 <- Pmax(x1, y1, z1)
pm2 <- pmax(x1, y1, z1, na.rm = TRUE)
all(pm1 == pm2)
benchr::benchmark(pmax(x1, y1, z1, na.rm = TRUE),
Pmax(x1, y1, z1))
Benchmark summary:
Time units : milliseconds
expr n.eval min lw.qu median mean up.qu max total relative
pmax(x1, y1, z1, na.rm = TRUE) 100 1.34 1.37 1.39 1.44 1.46 1.74 144 1.00
Pmax(x1, y1, z1) 100 13.30 13.50 13.80 19.90 15.70 67.50 1990 9.88
编辑：
我删除了一些循环，只是在 Rcpp 之外用 NA 替换了 -1，它加快了一点，但仍然没有超过 R base pmax。
虽然 Rcpp::pmax 是一个很好的实现，但它只处理两个向量，不确定它是否可以处理缺失值。当缺少值时，我得到了不同的结果。
第二次尝试是：
Pmax1 <- function(...) {
args_list <- list(...)
cppFunction("
NumericVector cpp_pmax(List args) {
List args0 = args[0];
int n_arg = args.length();
int n_vec = args0.length();
NumericVector out(n_vec);
out = args[0];
for (int i = 1; i < n_arg; ++i) {
NumericVector pa(n_vec);
pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (R_IsNA(out[j])) {
out[j] = -1;
}
if (R_IsNA(pa[j])) {
pa[j] = -1;
}
out[j] = std::max(out[j], pa[j]);
}
}
return out;
}
")
output <- cpp_pmax(args_list)
output[output == -1] <- NA
return(output)
}
Pmax2 <- function(...) {
args_list <- list(...)
cppFunction("
NumericVector cpp_pmax(List args) {
NumericVector out = args[0];
int n_arg = args.length();
int n_vec = out.length();
for (int j = 0; j < n_vec; ++j) {
if (NumericVector::is_na(out[j])) out[j] = -1;
}
for (int i = 1; i < n_arg; ++i) {
NumericVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (NumericVector::is_na(pa[j])) pa[j] = -1;
out[j] = std::max(out[j], pa[j]);
}
}
return out;
}
")
output <- cpp_pmax(args_list)
output[output == -1] <- NA
return(output)
}
n <- 200000
x <- sample(0:5, n, replace = TRUE)
y <- sample(0:5, n, replace = TRUE)
z <- sample(0:5, n, replace = TRUE)
w <- sample(0:5, n, replace = TRUE)
x[sample(1:n, 900)]<-NA
y[sample(1:n, 600)]<-NA
z[sample(1:n, 700)]<-NA
z[sample(1:n, 800)]<-NA
benchr::benchmark(pmax(x, y, z, w, na.rm = TRUE),
Pmax1(x, y, z, w),
Pmax2(x, y, z, w))
Benchmark summary:
Time units : milliseconds
expr n.eval min lw.qu median mean up.qu max total relative
pmax(x, y, z, w, na.rm = TRUE) 100 2.38 2.43 2.46 2.46 2.48 2.6 246 1.00
Pmax1(x, y, z, w) 100 16.00 16.90 17.20 19.40 17.70 61.2 1940 6.98
Pmax2(x, y, z, w) 100 9.44 9.74 9.90 11.30 10.10 45.6 1130 4.02
有没有人知道如何让它比 R base pmax 更快？
我们的想法是拥有一个通用函数来处理不同数量的向量，所有这些都在 Rcpp 函数中。
根据@DirkEddelbuettel 和@Cole 的回答更新
感谢您帮助优化代码。受@DirkEddelbuettel 和@Cole 回答的启发，我只是添加 Rcpp::pmax 来删除其中一个循环，它也有助于加快速度。
library("bench")
library("Rcpp")
cppFunction("
IntegerVector cpp_pmax1(List args) {
IntegerVector tmp = args[0];
IntegerVector out = clone(tmp);
int n_arg = args.length();
int n_vec = out.length();
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (pa[j] > out[j]) out[j] = pa[j];
}
}
return out;
}
")
cppFunction("
IntegerVector cpp_pmax2(List args) {
IntegerVector tmp = args[0];
IntegerVector out = clone(tmp);
int n_arg = args.length();
int n_vec = out.length();
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
out = pmax(out, pa);
}
return out;
}
")
Pmax1 <- function(...) {
cpp_pmax1(list(...))
}
Pmax2 <- function(...) {
cpp_pmax2(list(...))
}
n <- 200000
x <- sample(0:5, n, replace = TRUE)
y <- sample(0:5, n, replace = TRUE)
z <- sample(0:5, n, replace = TRUE)
w <- sample(0:5, n, replace = TRUE)
k <- sample(0:5, n, replace = TRUE)
x[sample(1:n, 900)] <- NA
y[sample(1:n, 600)] <- NA
z[sample(1:n, 700)] <- NA
w[sample(1:n, 800)] <- NA
k[sample(1:n, 800)] <- NA
pm0 <- pmax(x, y, z, w, k, na.rm = TRUE)
pm1 <- Pmax1(x, y, z, w, k)
pm2 <- Pmax2(x, y, z, w, k)
benchr::benchmark(pmax(x, y, z, w, k, na.rm = TRUE),
Pmax1(x, y, z, w, k),
Pmax2(x, y, z, w, k))
Benchmark summary:
Time units : microseconds
expr n.eval min lw.qu median mean up.qu max total relative
pmax(x, y, z, w, k, na.rm = TRUE) 100 2880 2900 2920 3050 3080 8870 305000 5.10
Pmax1(x, y, z, w, k) 100 2150 2180 2200 2310 2350 8060 231000 3.85
Pmax2(x, y, z, w, k) 100 527 558 572 812 719 7870 81200 1.00
谢谢！
顺便说一句，请注意 Rcpp 糖已经有 Rcpp::pmax():
> library(Rcpp)
> cppFunction("NumericVector pm(NumericVector x, NumericVector y) {
+ return pmax(x,y);}")
> pm(10.0*(1:10), rep(50, 10))
[1] 50 50 50 50 50 60 70 80 90 100
> pm(10.0*(1:10), c(rep(50, 8), NA, 50))
[1] 50 50 50 50 50 60 70 80 NA 100
>
可能还有另一个更通用的功能的范围，但希望这也可以作为基准对您有所帮助。
编辑： 在我的第一个版本中，当我打算调用 pm()（使用 Rcpp::pmax()）时，我不小心调用了 pmax()。结果一样。
pm() 和 pmax() 与人们预期的速度顺序大致相同，因为两者都是矢量化的：
> library(microbenchmark)
> set.seed(123)
> x <- cumsum(rnorm(1e6))
> y <- cumsum(rnorm(1e6))
> microbenchmark(pmax(x,y), pm(x,y))
Unit: milliseconds
expr min lq mean median uq max neval cld
pmax(x, y) 3.94342 4.07488 4.66378 4.15433 5.39961 7.81931 100 a
pm(x, y) 3.58781 3.68886 4.74249 3.75815 5.38444 22.31268 100 a
>
我想你可以尝试 fcoalesce + fifelse（都来自 data.table 包）来定义你的 Pmax 函数，如下所示
Pmax <- function(..., na.rm = FALSE) {
u <- list(...)
if (na.rm) {
return(
Reduce(function(x, y) {
x <- fcoalesce(x, y)
y <- fcoalesce(y, x)
fifelse(x <= y, y, x)
}, u)
)
}
Reduce(function(x, y) fifelse(x <= y, y, x), u)
}
Benchmark（使用 OP post 中的数据进行测试）
- 如果启用
na.rm = TRUE，Pmax比基本 Rpmax 稍慢
> microbenchmark::microbenchmark(
+ pmax(x1, y1, z1, na.rm = TRUE),
+ Pmax(x1, y1, z1, na.rm = TRUE),
+ check = "equivalent",
+ unit = "relati ..." ... [TRUNCATED]
Unit: relative
expr min lq mean median uq
pmax(x1, y1, z1, na.rm = TRUE) 1.000000 1.00000 1.000000 1.000000 1.000000
Pmax(x1, y1, z1, na.rm = TRUE) 1.428545 1.87539 1.974959 2.022579 2.094833
max neval
1.000000 100
1.387139 100
- 如果您使用默认的
na.rm选项，您会发现Pmax比基础 Rpmax 稍快
> microbenchmark::microbenchmark(
+ pmax(x1, y1, z1),
+ Pmax(x1, y1, z1),
+ check = "equivalent",
+ unit = "relative"
+ )
Unit: relative
expr min lq mean median uq max neval
pmax(x1, y1, z1) 1.387953 1.32482 1.053983 1.220124 1.143867 0.266205 100
Pmax(x1, y1, z1) 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 100
从 bench::mark 中可以看出内存分配似乎存在一些问题。
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2(x, y, z, w))
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc
## <bch:expr> <bch:t> <bch:t> <dbl> <bch:byt>
## 1 pmax(x, y, z, w, na.rm = TRUE) 5.79ms 6.28ms 157. 781.3KB
## 2 Pmax2(x, y, z, w) 39.56ms 54.48ms 19.7 9.18MB
内存强制
与基数 pmax() 相比，内存分配 是 10 倍。您的 rcpp 比较直接，因此这暗示存在某种强制。在查看示例数据时，您正在将整数向量发送到数字签名。这造成了代价高昂的胁迫。让我们更新签名和代码以期待 IntegerVectors。为此，我只是将所有内容从 NumericVector 更改为 IntegerVector。
expression min median `itr/sec` mem_alloc
<bch:expr> <bch:t> <bch:t> <dbl> <bch:byt>
1 pmax(x, y, z, w, na.rm = TRUE) 1.89ms 2.33ms 438. 781.3KB
2 Pmax2_int(x, y, z, w) 37.42ms 49.88ms 17.6 2.32MB
重新编译
OP 代码在较大的函数代码中包含 cppFunction。除非我们需要在每个循环中重新编译它，否则我们可以编译然后从 R 中调用编译后的代码。这是此数据集大小的最大性能提升。
cppFunction("
IntegerVector cpp_pmax_pre(List args) {
IntegerVector out = args[0];
int n_arg = args.length();
int n_vec = out.length();
for (int j = 0; j < n_vec; ++j) {
if (IntegerVector::is_na(out[j])) out[j] = -1;
}
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (IntegerVector::is_na(pa[j])) pa[j] = -1;
out[j] = std::max(out[j], pa[j]);
}
}
return out;
}
")
Pmax2_int_pre <- function(...) {
args_list <- list(...)
output <- cpp_pmax_pre(args_list)
output[output == -1] <- NA
return(output)
}
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2_int_pre(x, y, z, w))
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc
## <bch:expr> <bch:> <bch:> <dbl> <bch:byt>
## 1 pmax(x, y, z, w, na.rm = TRUE) 2.31ms 2.42ms 397. 781.3KB
## 2 Pmax2_int_pre(x, y, z, w) 2.48ms 3.55ms 270. 2.29MB
更多内存和小优化
最后，我们还有更多的内存分配。这暗示我们可以做更多 - 在这种情况下，我们应该在 rcpp 中更新 NA_REAL。相关的，我们可以优化一下循环赋值。
cppFunction("
IntegerVector cpp_pmax_final(List args) {
IntegerVector out = args[0];
int n_arg = args.length();
int n_vec = out.length();
for (int j = 0; j < n_vec; ++j) {
if (IntegerVector::is_na(out[j])) out[j] = -1;
}
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
// simplify logic; if the element is not na and is greater than the out, update out.
if (!IntegerVector::is_na(pa[j]) & pa[j] > out[j]) out[j] = pa[j];
}
}
// update now in Rcpp instead of allocating vectors in R
for (int i = 0; i < n_vec; i++) {
if(out[i] == -1) out[i] = NA_INTEGER;
}
return out;
}
")
Pmax2_final <- function(...) {
cpp_pmax_final(list(...))
}
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2_final(x, y, z, w))
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc
## <bch:expr> <bch:> <bch:> <dbl> <bch:byt>
## 1 pmax(x, y, z, w, na.rm = TRUE) 2ms 2.08ms 460. 781.3KB
## 2 Pmax2_final(x, y, z, w) 1.19ms 1.45ms 671. 2.49KB
我们做到了*！我确信可能会有一些小的优化 - 我们访问了 pa[j] 三次，因此可能值得分配给一个变量。
奖金 - NA_INTEGER
根据Rcpp for Everyone，NA_INTEGER 应该等于最小整数值-2147483648。使用这个，我们可以删除 NA 的替换，因为 我们可以在处理 int 数据类型时直接与 NA 进行比较。
在实现过程中，我还发现了前一部分的一个问题——我们需要克隆初始参数，这样我们就不会不小心通过引用更改它。不过，我们仍然比基础 pmax().
cppFunction("
IntegerVector cpp_pmax_last(List args) {
IntegerVector tmp = args[0];
IntegerVector out = clone(tmp);
int n_arg = args.length();
int n_vec = out.length();
for (int i = 1; i < n_arg; ++i) {
IntegerVector pa = args[i];
for (int j = 0; j < n_vec; ++j) {
if (pa[j] > out[j]) out[j] = pa[j];
}
}
return out;
}
")
Pmax2_last <- function(...) {
cpp_pmax_last(list(...))
}
bench::mark(pmax(x, y, z, w, na.rm = TRUE),
Pmax2_last(x, y, z, w),
)
## # A tibble: 2 x 13
## expression min median `itr/sec` mem_alloc `gc/sec`
## <bch:expr> <bch:> <bch:> <dbl> <bch:byt> <dbl>
## 1 pmax(x, y, z, w, na.rm = TRUE) 5.98ms 6.36ms 154. 781KB 0
## 2 Pmax2_last(x, y, z, w) 5.09ms 5.46ms 177. 784KB 0