关于摆脱循环的建议答案

【问题标题】：Advice wanted on getting rid of loops关于摆脱循环的建议
【发布时间】：2010-12-18 20:38:43
【问题描述】：

我编写了一个程序来处理 3n + 1 问题（又名“奇妙的数字”和其他各种问题）。但它有一个双循环。我怎么能矢量化它？

代码是

count <- vector("numeric", 100000)
L <- length(count)

for (i in 1:L)
{
x <- i
   while (x > 1)
   {
   if (round(x/2) == x/2) 
     {
     x <- x/2
     count[i] <- count[i] + 1 
     } else
     {
     x <- 3*x + 1
     count[i] <- count[i] + 1
     }
   }
}

谢谢！

【问题讨论】：

我将盗用这个作为 R 中令人尴尬的并行进程的示例！谢谢！

标签： r loops vectorization

【解决方案1】：

我通过创建一个向量 x 来实现这种“由内而外”，其中第 i 个元素是算法每次迭代后的值。结果比较容易理解为

f1 <- function(L) {
    x <- seq_len(L)
    count <- integer(L)
    while (any(i <- x > 1)) {
        count[i] <- count[i] + 1L
        x <- ifelse(round(x/2) == x/2, x / 2, 3 * x + 1) * i
    }
    count
}

这可以优化为 (a) 仅跟踪那些仍在运行的值（通过 idx）和 (b) 避免不必要的操作，例如，ifelse 对所有 x 值的两个参数求值，x/2 求值两次。

f2 <- function(L) {
    idx <- x <- seq_len(L)
    count <- integer(L)
    while (length(x)) {
        ix <- x > 1
        x <- x[ix]
        idx <- idx[ix]
        count[idx] <- count[idx] + 1L
        i <- as.logical(x %% 2)
        x[i] <- 3 * x[i] + 1
        i <- !i
        x[i] <- x[i] / 2
    }
    count
}

用 f0 原来的功能，我有

> L <- 10000
> system.time(ans0 <- f0(L))
   user  system elapsed 
  7.785   0.000   7.812 
> system.time(ans1 <- f1(L))
   user  system elapsed 
  1.738   0.000   1.741 
> identical(ans0, ans1)
[1] TRUE
> system.time(ans2 <- f2(L))
   user  system elapsed 
  0.301   0.000   0.301 
> identical(ans1, ans2)
[1] TRUE

一个调整是将奇数值更新为 3 * x[i] + 1，然后无条件地除以 2

x[i] <- 3 * x[i] + 1
count[idx[i]] <- count[idx[i]] + 1L
x <- x / 2
count[idx] <- count[idx] + 1

用这个作为 f3（不知道为什么今天早上 f2 变慢了！）我明白了

> system.time(ans2 <- f2(L))
   user  system elapsed 
   0.36    0.00    0.36 
> system.time(ans3 <- f3(L))
   user  system elapsed 
  0.201   0.003   0.206 
> identical(ans2, ans3)
[1] TRUE

似乎在除以二阶段可以采取更大的步骤，例如，8 是 2^3，所以我们可以采取 3 步（加 3 来计数）并完成，20 是 2^2 * 5所以我们可以分两步，在 5 处进入下一次迭代。实现？

【讨论】：

【解决方案2】：

因为您需要迭代 x 的值，所以您无法真正对其进行矢量化。在某些时候，R 必须分别依次处理 x 的每个值。您也许可以在不同的 CPU 内核上运行计算以加快速度，也许可以使用同名包中的foreach。

否则，（这只是对您隐藏循环），将循环的主体包装为一个函数，例如：

wonderous <- function(n) {
    count <- 0
    while(n > 1) {
        if(isTRUE(all.equal(n %% 2, 0))) {
            n <- n / 2
        } else {
            n <- (3*n) + 1
        }
        count <- count + 1
    }
    return(count)
}

然后您可以使用sapply() 对一组数字运行函数：

> sapply(1:50, wonderous)
 [1]   0   1   7   2   5   8  16   3  19   6  14   9   9  17  17
[16]   4  12  20  20   7   7  15  15  10  23  10 111  18  18  18
[31] 106   5  26  13  13  21  21  21  34   8 109   8  29  16  16
[46]  16 104  11  24  24

或者您可以使用Vectorize 返回wonderous 的矢量化版本，它本身就是一个函数，可以对您隐藏更多内容：

> wonderousV <- Vectorize(wonderous)
> wonderousV(1:50)
 [1]   0   1   7   2   5   8  16   3  19   6  14   9   9  17  17
[16]   4  12  20  20   7   7  15  15  10  23  10 111  18  18  18
[31] 106   5  26  13  13  21  21  21  34   8 109   8  29  16  16
[46]  16 104  11  24  24

~~我认为这是目前使用标准 R 工具所能达到的程度。~~@Martin Morgan 表明，通过巧妙地解决以下问题，您可以做得比这更好。确实使用了 R 的矢量化能力。

【讨论】：

我就是这样做的（但使用降雪包进行并行执行）。

【解决方案3】：

另一种方法可以识别一个人经常重新访问低数字，那么为什么不记住它们并节省重新计算的成本呢？

memo_f <- function() {
    e <- new.env(parent=emptyenv())
    e[["1"]] <- 0L
    f <- function(x) {
        k <- as.character(x)
        if (!exists(k, envir=e))
            e[[k]] <- 1L + if (x %% 2) f(3L * x + 1L) else f(x / 2L)
        e[[k]]
    }
    f
}

给了

> L <- 100
> vals <- seq_len(L)
> system.time({ f <- memo_f(); memo1 <- sapply(vals, f) })
   user  system elapsed 
  0.018   0.000   0.019 
> system.time(won <- sapply(vals, wonderous))
   user  system elapsed 
  0.921   0.005   0.930 
> all.equal(memo1, won) ## integer vs. numeric
[1] TRUE

这可能不会很好地并行化，但是对于 50 倍的加速，这可能不是必需的？此外，递归可能会变得太深，但可以将递归写成循环（无论如何这可能更快）。

【讨论】：