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 中实现这些方程，但这样做要具有合理的“数学可读性”水平。特别是，我不是在寻找最短或最快执行的代码块，而是直观地反映原始数学表达式的代码块。对于上面的示例，我知道计算向量 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 的理解，但对语法进行了某些修改。还发布了一些代码 here 和 here 以及一个仅限 github 的包 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)

更新

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