r - 除了使用 R 基函数之外，是否有一种有效的方法可以获得 "pmax"？

Question

我想使用 Rcpp 创建一个函数，它可以胜过 R base 中的 pmax 函数。
我还尝试处理 Rcpp 函数内部的缺失值，这可能不是一个好主意。
所有向量都必须有一些缺失值，并且它们都是正值。这就是我将缺失重新编码为 -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 基础 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
  

谢谢!

Answer 1

顺便说一句，请注意 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
> 

很可能还有另一个更通用的功能的空间，但希望这也可以帮助您作为基准。
编辑:在我的第一个版本中，我不小心调用了 pmax()当我打算打电话时 pm() (使用 Rcpp::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
>