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

问题描述

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

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 apprach 快得多）并且在数学上很直观；然而，代码读起来很笨拙（只有二维的七行）。因此，我想问一下：你能想出一个更好的替代方法来解决这个问题（以及更复杂的变体），同时又不牺牲太多的计算速度、数学直觉或代码可读性？

解决方法

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]) * R[i])))
## [1] -144 -240  -44 -144 -260 -112

4) 推导式 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)

更新

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

algebra r r