c++ - 如何加速 R 代码在特定方向上移动矩阵元素 - Rccp？

Question

我开发了代码，可以围绕大型矩阵移动元素以模拟水文条件。这部分代码是一个函数，在一个小时的模拟中被多次调用，而当我试图对 30 年的数据进行建模时，最终会耗费大量时间。

在名为 flowingrunoff 的原始矩阵(2710 行 x 7153 列)中，一些像元向东移动，一些向东南移动，一些向南移动等等(即到 8 个罗盘点之一)。这已经在早些时候确定并保存在 8 个相同大小的矩阵中，如果单元格朝那个方向移动则为 1，否则为 0(这些称为 fdrecells、fdrsecells、fdrscells 等)。前面还创建了一个由全 0 组成的矩阵，该矩阵比其他矩阵大两行两列(称为 emptymovematrix)。

现在转到移动单元格的函数，我称之为 route()。这是代码，它使用 set.submatrix 将矩阵移动到 emptymovematrix 的特定部分，如下所示:

route <- function() {    

flowingrunoff[is.na(flowingrunoff)] <- 0
flowingbaseflow[is.na(flowingbaseflow)] <- 0

fdrerunoff <- fdrecells*flowingrunoff
movefdrerunoff <- set.submatrix(emptymovematrix,fdrerunoff,2,3)

fdrserunoff <- fdrsecells*flowingrunoff
movefdrserunoff <- set.submatrix(emptymovematrix,fdrserunoff,3,3)

fdrsrunoff <- fdrscells*flowingrunoff
movefdrsrunoff <- set.submatrix(emptymovematrix,fdrsrunoff,3,2)

fdrswrunoff <- fdrswcells*flowingrunoff
movefdrswrunoff <- set.submatrix(emptymovematrix,fdrswrunoff,3,1)

fdrwrunoff <- fdrwcells*flowingrunoff
movefdrwrunoff <- set.submatrix(emptymovematrix,fdrwrunoff,2,1)

fdrnwrunoff <- fdrnwcells*flowingrunoff
movefdrnwrunoff <- set.submatrix(emptymovematrix,fdrnwrunoff,1,1)

fdrnrunoff <- fdrncells*flowingrunoff
movefdrnrunoff <- set.submatrix(emptymovematrix,fdrnrunoff,1,2)

fdrnerunoff <- fdrnecells*flowingrunoff
movefdrnerunoff <- set.submatrix(emptymovematrix,fdrnerunoff,1,3)

因此创建了 8 个新矩阵。现在可以将这些组合起来，以求出单元移动后的径流总和，然后将其修剪为原始矩阵的大小。

movedrunoff <- Reduce(function(x,y){x + y}, list(movefdrerunoff,movefdrserunoff,movefdrsrunoff,movefdrswrunoff,movefdrwrunoff,movefdrnwrunoff,movefdrnrunoff,movefdrnerunoff))

movedrunoff <- movedrunoff[2:(nrow(movedrunoff)-1),2:(ncol(movedrunoff)-1)]

}

这有效并且比我原来使用 rbind 和 cbind 的方法快得多(3 到 4 倍)，如下所示:

fdrerunoff <- fdrecells*flowingrunoff
movefdrerunoff <- cbind(0,0,fdrerunoff)
movefdrerunoff <- rbind(0,movefdrerunoff,0)

然而，尽管这段代码(连同一些围绕基本流程移动的类似代码)每个周期只需要 6 到 10 秒，但每小时模拟了 5 到 10 个周期——你可以看到这在时间上加起来很长.整个代码每小时大约需要 70 秒，而 30 多年需要大约 212 天才能完成。有关信息，route() 函数是从另一个名为 routing() 的函数调用的，如下所示:

routing <- function () {

final <- replicate(replication,route())

}

那么，关于如何改进此代码或策略的任何想法。我试过 shift.right(在 Matrix 中)，我试过转换为栅格并移动范围，我试过 rbind 和 cbind 但都比上述方法慢得多。我相信加快速度的真正关键是在 C++ 中使用调用代码，可能使用 Rccp，但不知道从哪里开始。所以有任何更巧妙的 R 技术可以移动矩阵，或者如果您知道如何用 C++ 编写代码，请告诉我。

令我惊讶的是，R 中的矩阵乘法非常快，但移动或复制它们或使用矩阵 [] 之类的东西相对较慢。请注意，我使用的是安装了英特尔数学核心函数库的 Revolutionary R。

非常感谢

安东尼

Answer 1

我建议使用邻接矩阵/转换矩阵方法。这样你就可以从 %*% 中获取能量

创建两个转换矩阵 transition1 <- diag(7153) transition2 <- diag(2710)

第一个控制一个轴上的转换，第二个控制另一个轴上的转换

然后为每个转换编辑要进行的转换的行/列 your_matrix %*% transition1和 transition2 %*% your_matrix应该给你移动的值。

编辑:如果你这样做，顺序很重要，如果你想同时做，那么你需要捕获reshape2，将你的原始data.frame融化成长格式，并制作一个方阵(2710*7153) 这可能对您的需求来说太大了。但是，如果您要对流进行建模，则可以避免使用这种两步近似

EDIT2:如果你只是做transition2 %*% your_matrix + your_matrix %*% transition1，你实际上可以绕过订单只要你确保你没有失去任何流量

EDIT3:所以，这是一个 6x8 的小例子(我假设你已经把它分解成你的 8 个方向，我假设你不允许流体留在原处)

basemat <- matrix(rep(c(1,1,0,0,0,0,0,0),6), nrow = 6)

transition_east <- diag(8)
transition_east <- transition_east[c(2:8,1),]

所以这里的转换矩阵被设置为做一些非常简单的事情，它将所有列向东移动一个。 (它包裹在边缘)

然后尝试: basemat

basemat %*% transition_east

basemat %*% transition_east %*% transition_east

等...