为什么一些逻辑运算符如此缓慢？

Question

在尝试优化我的代码时，我发现与 integer 或 numeric 上类似的操作相比，一些 logical 操作速度较慢。

所以我尝试将基本的布尔运算符 !、&、|、xor 重写如下：

my.not <- function(x) as.logical(1L - as.integer(x))
my.and <- function(e1, e2) as.logical(as.integer(e1) * as.integer(e2))
my.or  <- function(e1, e2) as.logical(as.integer(e1) + as.integer(e2))
my.xor <- function(e1, e2) as.logical(as.integer(e1) + as.integer(e2) == 1L)

测试一切是否按预期工作：

a <- sample(c(TRUE, FALSE), 1e6, TRUE)
b <- sample(c(TRUE, FALSE), 1e6, TRUE)

identical(!a, my.not(a))             # TRUE
identical(a & b, my.and(a, b))       # TRUE
identical(a | b, my.or(a, b))        # TRUE
identical(xor(a, b), my.xor(a, b))   # TRUE

现在正在进行基准测试:

library(microbenchmark)
microbenchmark(!a, my.not(a),
               a & b, my.and(a, b),
               a | b, my.or(a, b),
               xor(a, b), my.xor(a, b))

# Unit: milliseconds
#          expr       min        lq    median         uq       max neval
#            !a  1.237437  1.459042  1.463259   1.492671  17.28209   100
#     my.not(a)  6.018455  6.263176  6.414515  15.291194  70.16313   100
#         a & b 32.318530 32.667525 32.769014  32.973878  50.55528   100
#  my.and(a, b)  8.010022  8.592776  8.750786  18.145590  78.38736   100
#         a | b 32.030545 32.383769 32.506937  32.820720 102.43609   100
#   my.or(a, b) 12.089538 12.434793 12.663695  22.046841  32.19095   100
#     xor(a, b) 94.892791 95.480200 96.072202 106.104000 164.19937   100
#  my.xor(a, b) 13.337110 13.708025 14.048350  24.485478  29.75883   100

从结果来看，! 运算符似乎是唯一一个能够与我的实现相当的。而另外三个运算符则慢了几倍，这对于 Primitive 函数来说有点尴尬。我甚至希望一个良好实现的布尔运算符应该比整数运算（我的函数实现方式）要快得多。

问题：为什么？实现不好吗？还是原始函数做了一些好事情（例如错误检查、特殊情况），而我的函数没有做到？

- flodel

---

不要忘记NA，例如NA|TRUE、NA|FALSE和名称c(x=TRUE)。 - Martin Morgan

Answer 1

稍微看一下C语言的实现，逻辑和数学操作以不同的方式实现它们的循环。逻辑操作会执行类似以下的操作(在logic.c:327)

library(inline)

or1 <- cfunction(c(x="logical", y="logical"), "
    int nx = LENGTH(x), ny = LENGTH(y), n = nx > ny ? nx : ny;
    SEXP ans = PROTECT(allocVector(LGLSXP, n));
    int x1, y1;
    for (int i = 0; i < n; i++) {
        x1 = LOGICAL(x)[i % nx];
        y1 = LOGICAL(y)[i % ny];
        if ((x1 != NA_LOGICAL && x1) || (y1 != NA_LOGICAL && y1))
            LOGICAL(ans)[i] = 1;
        else if (x1 == 0 && y1 == 0)
            LOGICAL(ans)[i] = 0;
        else
            LOGICAL(ans)[i] = NA_LOGICAL;
    }
    UNPROTECT(1);
    return ans;
")

每次迭代中都有两个模运算符%。相比之下，算术运算（在Itermacros.h:54中）执行类似于以下操作：

or2 <- cfunction(c(x="logical", y="logical"), "
    int nx = LENGTH(x), ny = LENGTH(y), n = nx > ny ? nx : ny;
    SEXP ans = PROTECT(allocVector(LGLSXP, n));
    int x1, y1, ix=0, iy=0;
    for (int i = 0; i < n; i++) {
        x1 = LOGICAL(x)[ix];
        y1 = LOGICAL(x)[iy];
        if (x1 == 0 || y1 == 0)
            LOGICAL(ans)[i] = 0;
        else if (x1 == NA_LOGICAL || y1 == NA_LOGICAL)
            LOGICAL(ans)[i] = NA_LOGICAL;
        else
            LOGICAL(ans)[i] = 1;

        if (++ix == nx) ix = 0;
        if (++iy == ny) iy = 0;
    }
    UNPROTECT(1);
    return ans;
")

进行两个身份验证测试。这是一种跳过NA测试的版本。

or3 <- cfunction(c(x="logical", y="logical"), "
    int nx = LENGTH(x), ny = LENGTH(y), n = nx > ny ? nx : ny;
    SEXP ans = PROTECT(allocVector(LGLSXP, n));
    int x1, y1, ix=0, iy=0;
    for (int i = 0; i < n; ++i) {
        x1 = LOGICAL(x)[ix];
        y1 = LOGICAL(y)[iy];
        LOGICAL(ans)[i] = (x1 || y1);
        if (++ix == nx) ix = 0;
        if (++iy == ny) iy = 0;
    }
    UNPROTECT(1);
    return ans;
")

然后提供一个避免使用LOGICAL宏的版本

or4 <- cfunction(c(x="logical", y="logical"), "
    int nx = LENGTH(x), ny = LENGTH(y), n = nx > ny ? nx : ny;
    SEXP ans = PROTECT(allocVector(LGLSXP, n));
    int *xp = LOGICAL(x), *yp = LOGICAL(y), *ansp = LOGICAL(ans);
    for (int i = 0, ix = 0, iy = 0; i < n; ++i)
    {
        *ansp++ = xp[ix] || yp[iy];
        ix = (++ix == nx) ? 0 : ix;
        iy = (++iy == ny) ? 0 : iy;
    }
    UNPROTECT(1);
    return ans;
")

以下是一些时间表

microbenchmark(my.or(a, b), a|b, or1(a, b), or2(a, b), or3(a, b), or4(a, b))
Unit: milliseconds
        expr       min        lq    median       uq      max neval
 my.or(a, b)  8.002435  8.100143 10.082254 11.56076 12.05393   100
       a | b 23.194829 23.404483 23.860382 24.30020 24.96712   100
   or1(a, b) 17.323696 17.659705 18.069139 18.42815 19.57483   100
   or2(a, b) 13.040063 13.197042 13.692152 14.09390 14.59378   100
   or3(a, b)  9.982705 10.037387 10.578464 10.96945 11.48969   100
   or4(a, b)  5.544096  5.592754  6.106694  6.30091  6.94995   100

a|b和or1的区别反映了这里未实现的内容，比如属性、尺寸和对象的特殊处理。从or1到or2反映了不同回收方式的成本；我对这里存在差异感到惊讶。从or2到or3是NA安全性的成本。很难确定or4中额外的加速是否会在基础R实现中看到 - 在用户C代码中，LOGICAL()是一个宏，但在基础R中它是一个内联函数调用。

该代码使用-O2标志进行编译。

> system("clang++ --version")
Ubuntu clang version 3.0-6ubuntu3 (tags/RELEASE_30/final) (based on LLVM 3.0)
Target: x86_64-pc-linux-gnu
Thread model: posix

my.or 的运行时间在不同的独立R会话之间并不一致，有时会花费更长的时间；我不确定为什么会这样。上面的时间是使用 R version 2.15.3 Patched (2013-03-13 r62579) 运行时所用的时间；当前的 R-devel 版本看起来快了约10%。

- Martin Morgan

---

感谢你深入挖掘R语言内部的历史，这总是很有启发性的，可以让我们更好地了解R语言的内部实现。 - Kevin Ushey

也许值得在r-devel上提出。 - hadley

Answer 2

虽然我非常喜欢你的方法并且喜欢速度的提升，但是很遗憾当e1和e2的结构比向量更加复杂时，它们就失去了效果。

> dim(a) <- c(1e2, 1e4)
> dim(b) <- c(1e2, 1e4)
> 
> identical(!a, my.not(a))            
[1] FALSE
> identical(a & b, my.and(a, b))       
[1] FALSE
> identical(a | b, my.or(a, b))        
[1] FALSE
> identical(xor(a, b), my.xor(a, b))   
[1] FALSE

结构/维度

逻辑功能保留了结构和属性，这是很昂贵的，但也是有价值的。

T <- TRUE; F <- FALSE
A <- matrix(c(T, F, T, F), ncol=2)
B <- matrix(c(T, F, F, T), ncol=2)

> A & B
      [,1]  [,2]
[1,]  TRUE FALSE
[2,] FALSE FALSE

> my.and(A, B)
[1]  TRUE FALSE FALSE FALSE

缺失值处理

正如评论中指出的那样，还需要考虑缺失值（NA）-- 即需要更多的开销。

a <- c(T, F, NA, T)
b <- c(F, NA, T, T)
>     identical(!a, my.not(a))    
[1] TRUE
>     identical(a & b, my.and(a, b))       
[1] FALSE
>     identical(a | b, my.or(a, b))        
[1] FALSE
>     identical(xor(a, b), my.xor(a, b))   
[1] TRUE

属性

a <- c(T, F, NA, T)
b <- c(F, NA, T, T)

names(a) <- names(b) <- LETTERS[23:26]

> a & b
    W     X     Y     Z 
FALSE FALSE    NA  TRUE 

> my.and(a, b)
[1] FALSE    NA    NA  TRUE

---

不是速度问题

当然，说到这些，你的函数确实提供了很大的增益！如果你知道不必担心NAs等问题，也不在意结构，那么为什么不直接使用它们呢！

- Ricardo Saporta

Answer 3

此外，不要忘记在R中，逻辑值被存储为整数在内存中。因此，速度的提升可能并不会很大。

> a=rep(TRUE,1000)
> b=rep(1L,1000)
> c=rep(1,1000)
> class(a)
[1] "logical"
> class(b)
[1] "integer"
> class(c)
[1] "numeric"
> object.size(a)
4040 bytes
> object.size(b)
4040 bytes
> object.size(c)
8040 bytes

- Tiago Zortea

---

尝试执行 rep(TRUE, 1000), rep(1L, 1000), rep(1, 1000)。 - flodel

是的，我说了整数，但在我的代码中，我进行了数字比较，我扩展了代码以澄清。 - Tiago Zortea