R 中的代数：矢量化与 foreach 与 sapply

Question

我有一系列包含向量、矩阵和数组的和与乘积的方程组，例如这个：

Y_i = \sum_{s=1}^S (1-alpha_{i,s})*R_i,

其中 Y 和 R 是长度为 I 且元素分别为 Y_i 和 R_i 的向量，而 alpha 是一个具有 I 行和 S 列的矩阵。

现在我想在 R 中实现这些方程式，但要以 'mathematical readability' 的合理水平实现。特别是，我不是在寻找最短或执行速度最快的代码块，而是寻找能够直观地反映原始数学表达式的代码块。对于上面的示例，我知道计算向量 Y 的一种快速简便的方法是向量化：

Y <- rowSums((1-alpha)*R)

然而，考虑到更复杂的表达式、更多的操作和更多的维度，我发现使用 foreach 遍历所涉及的维度的循环来基本上复制纸上的方程更直观：

library(foreach)
Y <- foreach(i = 1:I, .combine = c) %:%
    foreach(s = 1:S, .combine = sum) %do% {
        (1-alpha[i,s])*R[i]
    }

我很喜欢这里的结构和.combine 参数，代码还是有些简洁的。不幸的是，这种方法的性能很糟糕，令人遗憾的是它不可行。然后我尝试了 sapply 循环：

Y <- sapply(1:I, function(i) {
    sum(
        sapply(1:S, function(s) {
            (1-alpha[i,s])*R[i]
        })
    )
})

这种方法既快速（不如矢量化方法快，但比 foreach 方法快得多）并且在数学上很直观；然而，代码阅读起来非常笨拙（七行代码只有两个维度）。因此我想问：你能想出一种更好的替代方法来解决这个问题（及其更复杂的变体）而不牺牲太多的计算速度、数学直觉或代码可读性吗？

Answer 1

1) for 仅对内部循环进行向量化会得到非常接近原始循环的结果。（我们使用最后注释中的输入。）

I <- nrow(alpha)
Y <- numeric(I)
for(i in 1:I) Y[i] <- sum((1 - alpha[i, ]) * R[i])
## [1] -144 -240  -44 -144 -260 -112

2) sapply 这也适用于类似的方法：

I <- nrow(alpha)
Y <- sapply(1:I, function(i) sum((1 - alpha[i, ]) * R[i]))
## [1] -144 -240  -44 -144 -260 -112

3) fn$ 使用 gsubfn 包中的 fn$ 作为函数的前缀将允许传入参数的函数指定为公式，因此我们可以这样写：

library(gsubfn)

I <- nrow(alpha)
S <- ncol(alpha)
fn$sapply(1:I, i ~ sum(fn$sapply(1:S, s ~ (1 - alpha[i, s]) * R[i])))
## [1] -144 -240  -44 -144 -260 -112

或者为了更简洁，我们定义 iter 如图所示，并且还使用 gsubfn 的多重赋值工具一次定义 I 和 S。

library(gsubfn)

iter <- fn$sapply
list[I, S] <- dim(alpha)

iter(1:I, i ~ sum(iter(1:S, s ~ (1 - alpha[i, s]) * R[i])))
## [1] -144 -240  -44 -144 -260 -112

4) Comprehensions CRAN 上有 3 个包支持 python-like comprehensions，并对语法进行了某些修改。还有一些代码已发布 here and here and a github only package lc here，我们不会对其进行审查。这些软件包按字母顺序列在下面。

4a) 理解

library(comprehenr)
packageVersion("comprehenr") # be sure to use version 0.6.9 or later

I <- nrow(alpha)
to_vec(for(i in 1:I) sum((1-alpha[i, ])*R[i]))
## [1] -144 -240  -44 -144 -260 -112

或者有两个索引：

I <- nrow(alpha)
J <- ncol(alpha)
to_vec(for(i in 1:I) sum(to_vec(for(j in 1:J) (1-alpha[i, j])*R[i])))
## [1] -144 -240  -44 -144 -260 -112

4b) eList 有一个新的包，eList，支持列表和向量理解。这个包的一个显着特点（这里没有显示）是它支持并行处理，参数只有很小的变化。

library(eList)
packageVersion("eList") # be sure to use version 0.2,0 or later

I <- nrow(alpha)
Num(for (i in 1:I) sum((1-alpha[i, ])*R[i]))
## [1] -144 -240  -44 -144 -260 -112

或同时使用 i 和 s:

library(eList)

I <- nrow(alpha)
S <- ncol(alpha)

Num(for(i in 1:I) Sum(for(s in 1:S) (1-alpha[i,s]) * R[i]))
## [1] -144 -240  -44 -144 -260 -112

4c) listcompr 这是另一个支持理解的包。它的语法与上述两个略有不同，更接近 Python.

library(listcompr)

I <- nrow(alpha)
gen.vector(sum((1-alpha[i, ])*R[i]), i = 1:I)
## [1] -144 -240  -44 -144 -260 -112

或同时使用两个索引：

I <- nrow(alpha)
J <- ncol(alpha)
gen.vector(sum(gen.vector((1-alpha[i, j])*R[i], j = 1:J)), i = 1:I)
## [1] -144 -240  -44 -144 -260 -112

5) nimble 如果以上方法不够快，我们可以考虑使用 nimble 包，它将类似 R 的代码和一些类型定义翻译成 C++。

library(nimble)

calc <- nimbleFunction(
  run = function(alpha = double(2), R = double(1)) {
    I <- dim(alpha)[1]
    Y <- numeric(I)
    for(i in 1:I) Y[i] <- sum((1 - alpha[i, ]) * R[i])
    return(Y)
    returnType(double(1))
  }
)

Ccalc <- compileNimble(calc)

# test
Ccalc(alpha, R)
## [1] -144 -240  -44 -144 -260 -112

6) einsum einsum 包支持爱因斯坦张量符号。第一个参数的左侧用逗号分隔成两组，每组在后续参数的一个输入中定义索引。右侧的索引是输出的相应索引。此包具有生成 C++ 代码然后执行它的能力（此处未显示）。

library(einsum)

einsum("ij,i -> i", 1-alpha, R)
## [1] -144 -240  -44 -144 -260 -112

备注

用于测试的一些输入：

alpha <- matrix(1:24, 6)
R <- c(4, 6, 1, 3, 5, 2)

更新

根据新版本和其他发现重新安排了演示文稿，添加了其他方法并更新了理解部分。

R 中的代数：矢量化与 foreach 与 sapply

Algebra in R: Vectorization vs. foreach vs. sapply

r

algebra

备注

更新