非常感谢 @Mikko Marttila 更正了函数 3 和 4 并提供了函数 5 的想法。
最好使用矢量化选项而不是显式循环来处理 R。比如k的内循环：
for (k in 1:Vcut) {
b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
}
这和说的一样
(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V)
这个小改动让我们这部分函数的性能提升了 500 倍：
Unit: microseconds
expr min lq mean median uq max neval
k_loop 13186.7 13603.2 14605.471 13832.9 14517.8 41935.1 100
k_vectorized 16.4 17.6 25.559 28.8 32.0 52.7 100
现在，如果我们用 i 查看外部循环，我们会发现实际上没有必要逐行循环。我们可以为 sum(b[k]) 语句创建一个矩阵，将其变为：
(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V)
进入这个：
(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V)
这只是为我们节省了N*N*k 循环。在您的情况下，这是 646,400 个循环。
总而言之，我们会：
Arbfunc3 <- function(dV){
for (n in 1:Vcut) {
sum_b = colSums((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V))
fNC[, n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
fNDC[, n] = DC[,n]/fNC[,n]
}
}
对于这个替代方案，我进行微基准测试的中位时间是 750 毫秒。
为了进一步提高性能，我们需要解决V[n] - V。值得庆幸的是，R 有一个函数 - outer(V, V, '-')，这将生成一个包含我们需要的所有组合的矩阵。
Arbfunc4 <- function(dV) {
sum_b = apply((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(outer(V, V, '-')) / dV) / V, 2, function(x) colSums(x * t(C)))
fNC = exp(1*abs(V))*(1/(2*dV))*t(sum_b)
fNDC= DC/t(fNC)
fNDC
}
感谢@Mikko Marttila 提出的摆脱应用点积的建议。
Arbfunc5 <- function(dV) {
a = (V[2] - V[1]) * exp(-abs(V)) * t(C) / V
b = exp(abs(outer(V, V, "-")) / dV) %*% a
fNC = exp(1*abs(V))*(1/(2*dV))*(b)
fNDC= DC/t(fNC)
fNDC
}
这是每个解决方案的 system.time（Arbfunc2 是消除 k_loop）。优化后的解决方案比原始解决方案快 2,600 倍。
> system.time(Arbfunc(0.5))
user system elapsed
78.03 0.39 79.72
> system.time(Arbfunc2(0.5))
user system elapsed
10.41 0.03 10.46
> system.time(Arbfunc3(0.5))
user system elapsed
0.69 0.13 0.81
> system.time(Arbfunc4(0.5))
user system elapsed
0.43 0.05 0.47
> system.time(Arbfunc5(0.5))
user system elapsed
0.03 0.00 0.03
最终编辑：这是我在重新启动 R 并清空环境后运行的完整代码。没有错误：
## subsitutes for original data
DC <- matrix(rnorm(10), ncol=101, nrow=6400)
C <- matrix(rnorm(20), ncol=101, nrow=6400)
N <- 80
Vcut <- ncol(DC)
V <- seq(-2.9,2.5,length=Vcut)
# Unneeded for Arbfunc4 adn Arbfunc5
# Corrected from NA to NA_real_ to prevent coercion from logical to numeric
# h/t to @HenrikB
fNC <- matrix(NA_real_, nrow=(N*N), ncol=Vcut)
fNDC <- matrix(NA_real_, nrow=(N*N), ncol=Vcut)
Arbfunc <- function(dV){
b <- matrix(NA, nrow=1, ncol=Vcut)
for(i in 1:(N*N)) {
for (n in 1:Vcut) {
for (k in 1:Vcut) {
b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
}
fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum(b[]))
fNDC[i,n] = DC[i,n]/fNC[i,n]
}
}
fNDC
}
Arbfunc2 <- function(dV){
b <- matrix(NA, nrow=1, ncol=Vcut)
for(i in 1:(N*N)) {
for (n in 1:Vcut) {
sum_b = sum((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V))
fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
fNDC[i,n] = DC[i,n]/fNC[i,n]
}
}
fNDC
}
Arbfunc3 <- function(dV){
for (n in 1:Vcut) {
sum_b = colSums((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V))
fNC[, n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
fNDC[, n] = DC[,n]/fNC[,n]
}
fNDC
}
Arbfunc4 <- function(dV) {
sum_b = apply((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(outer(V, V, '-')) / dV) / V, 2, function(x) colSums(x * t(C)))
fNC = exp(1*abs(V))*(1/(2*dV))*t(sum_b)
DC/t(fNC)
}
Arbfunc5 <- function(dV) {
#h/t to Mikko Marttila for dot product
a = (V[2] - V[1]) * exp(-abs(V)) * t(C) / V
b = exp(abs(outer(V, V, "-")) / dV) %*% a
fNC = exp(1*abs(V))*(1/(2*dV))*(b)
DC/t(fNC)
}
#system.time(res <- Arbfunc(0.5))
system.time(res2 <- Arbfunc2(0.5))
system.time(res3 <- Arbfunc3(0.5))
system.time(res4 <- Arbfunc4(0.5))
system.time(res5 <- Arbfunc5(0.5))
all.equal(res2,res3,res4,res5)
正如@HenrikB 提到的，fNC 和fNDC 初始化为逻辑矩阵。这意味着我们在将它们强制为real 矩阵时会受到性能影响。对这个数据集执行不正确的操作是 1 毫秒的一次性命中，但如果这种强制处于循环中，它真的可以加起来。
mat_NA_real_ <- function() {
mat = matrix(NA_real_, nrow = 6400, ncol = 101)
mat[1,1] = 1
}
mat_NA <- function() {
mat = matrix(NA, nrow = 6400, ncol = 101)
mat[1,1] = 1
}
microbenchmark(mat_NA_real_(), mat_NA())
Unit: microseconds
expr min lq mean median uq max neval
mat_NA_real_() 979.5 992.25 1490.081 998.65 1021.1 7612.5 100
mat_NA() 1865.8 1883.30 3793.119 1911.30 5335.4 53635.2 100