我有一系列方程组，其中包含向量、矩阵和数组的和和乘积，例如：
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]
})
)
})
这种方法既快(不如矢量化方法快，但比foreachapprach快得多)和数学上直观;然而，代码阅读起来非常笨拙(只有两个维度的七行)。因此，我想问：你能想到一种更好的替代方法来解决这个问题(及其更复杂的变体)而不牺牲太多的计算速度、数学直觉或代码可读性吗？
1)对于矢量化，仅内部循环将给出非常接近原始循环的内容。 (我们使用最后注释中的输入。
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)理解CRAN 上有 3 个包支持类似 python 的推导，并对语法进行了某些修改。 这里和这里也发布了一些代码，这里还有一个只有 github 的包 lc，我们不会审查。 下面按字母顺序列出了这些软件包。
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)灵活如果上述方法不够快，我们可以考虑使用灵活的包，它将类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 包支持 Einstein张量符号。 第一个参数的左侧用逗号分隔为两组，每组定义后续参数的一个输入中的索引。 右侧的索引是输出的相应索引。 此包能够生成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)
更新
根据新版本和额外的发现，重新安排了演示文稿，添加了其他方法并更新了有关理解的部分。