6

我正在优化我的代码，我遇到了一些问题。我知道 R 中最大的加速来自矢量化代码而不是使用循环。但是，我的数据在列表中，我不确定是否可以矢量化我的代码。我尝试过使用这些apply功能（例如lapply，vapply），但我读到这些函数只是为了编写更简洁的代码，实际上是在底层使用循环！

这是我的代码中的三个最大瓶颈，尽管我认为第一部分没有什么可以做的。

1) 读取数据

我使用 1000 个尺寸为 277x349 的矩阵批量工作。doMC这是我脚本中最大的瓶颈，但我通过使用该包来利用该foreach功能的多个内核，稍微缓解了这个问题。这会产生一个包含 1000 个 277x349 矩阵的列表。

出于问题的目的，假设我们有 1000 个尺寸为 277 x 349 的矩阵的列表

# Fake data
l <- list()
for(i in 1:1000) {
  l[[i]] <- matrix(rnorm(277*349), nrow=277, ncol=349)
}

2) 瓶颈 #1

我需要与一些参考矩阵（相同尺寸）进行比较。这导致将列表中的 1000 个矩阵与我的参考矩阵进行比较，以获得 1000 个距离的向量。如果我知道矩阵具有相同的维度，我可以向量化这一步吗？

这是一些代码：

# The reference matrix
r <- matrix(rnorm(277*349), nrow=277, ncol=349)
# The number of non NA values in matrix. Do not need to worry about this...
K <- 277*349

# Make a function to calculate distances
distance <- function(xi, xj, K, na.rm=TRUE) {
  sqrt(sum((xi - xj)^2, na.rm=na.rm)/K)
}

# Get a vector containing all the distances
d <- vapply(l, distance, c(0), xj=r, K=K)

这一步是可以忍受的快速使用vapply，但它是代码中第三慢的部分。

3) 瓶颈 #2

我现在想为我的参考矩阵制作 J 个“最接近”矩阵的加权平均矩阵。（有一个排序步骤，但为简单起见假设d[1] < d[2] < ... < d[1000]）。我想获得 J=1,2,...,1000 时的加权平均矩阵

# Get the weighted matrix
weightedMatrix <- function(listOfData, distances, J) {
  # Calculate weights:
  w <- d[1:J]^{-2} / sum(d[1:J]^{-2})

  # Get the weighted average matrix
  # *** I use a loop here ***
  x_bar <- matrix(0, nrow=nrow(listOfData[[1]]), ncol=ncol(listOfData[[1]]))
  for(i in 1:J) {
    x_bar <- x_bar + {listOfData[[i]] * w[i]}
  }

  return(x_bar)
}

# Oh no! Another loop...
res <- list()
for(i in 1:length(l) ) {
  res[[i]] <- weightedMatrix(l, d, J=i)
}

我有点难过。我没有看到一种直观的方法来对矩阵列表进行矢量化操作。

我正在编写的脚本会被相当频繁地调用，所以即使是一点点改进也可以加起来！

编辑：

RE: 1) 读取数据

我忘了说我的数据是特殊格式的，所以我必须使用特殊的数据读取函数来读取R中的数据。文件是netcdf4格式，我正在使用nc_open包中的函数ncdf4来访问文件，然后我必须使用该ncvar_get函数来读取感兴趣的变量。好处是可以从磁盘读取文件中的数据，然后我可以将数据读入内存，ncvar_get然后用 R 对它们进行操作。

话虽如此，虽然我知道我的矩阵的大小以及我将拥有多少个矩阵，但我还是用数据列表提出了我的问题，因为foreach使我能够进行并行计算的函数输出了并行化循环的结果一个列表。我发现使用该foreach功能，数据读取步骤大约快 3 倍。

我想我可以在之后将数据排列为 3d 数组，但也许分配 3d 数组所花费的时间可能比它节省的时间更多？明天我得试试。

编辑2：

以下是我对脚本的一些时间安排。

原始脚本：

[1] "Reading data to memory"
user  system elapsed 
176.063  44.070  26.611 

[1] "Calculating Distances"
user  system elapsed 
2.312   0.000   2.308 

[1] "Calculating the best 333 weighted matrices"
user  system elapsed 
63.697  28.495   9.092

到目前为止，我做了以下改进：（1）在读取数据之前预先分配列表，（2）根据 Martin Morgan 的建议改进了加权矩阵计算。

[1] "Reading data to memory"
user  system elapsed 
192.448  38.578  27.872 

[1] "Calculating Distances"
user  system elapsed 
2.324   0.000   2.326 

[1] "Calculating all 1000 weighted matrices"
user  system elapsed 
1.376   0.000   1.374

一些注意事项：

我在foreach循环中使用 12 个内核来读取数据（registerDoMC(12) )。整个脚本在改进之前/之后运行大约需要 40 秒 / 36 秒。

我的瓶颈 #2 的时间已经改进了很多。以前，我只计算前三分之一（即 333）的加权矩阵，但现在脚本可以在原始时间的一小部分内计算所有加权矩阵。

感谢您的帮助，稍后我将尝试调整我的代码，看看是否可以更改我的脚本以使用 3D 数组而不是列表。我现在要花一些时间来验证计算，以确保它们有效！

4

1 回答 1

13

我的“低悬的果实”（scan; 预分配和填充）似乎无关紧要，所以......

距离计算中的操作对我来说看起来已经足够矢量化了。可能您可以通过对所有矩阵进行一次距离计算来挤出一些额外的速度，但这可能会使代码更难理解。

weightedMatrix 计算看起来还有改进的余地。让我们计算一下

w <- d^(-2) / cumsum(d^(-2))

对于加权矩阵m，我认为连续矩阵之间的关系只是m' = m * (1 - w[i]) + l[[i]] * w[i]，所以

res <- vector("list", length(l))
for (i in seq_along(l))
    if (i == 1L) {
        res[[i]] = l[[i]] * w[[i]]
    } else  {
        res[[i]] = res[[i - 1]] * (1 - w[[i]])  + l[[i]] * w[[i]]
    }

这将计算res从二次变为线性。我关于优于线性性能的想法只是一种（可能也是被误导的）预感。我没有追求过。

回到预分配和填充以及@flodel 的评论，我们有

f0 <- function(n) {
    ## good: pre-allocate and fill
    l = vector("list", n)
    for (i in seq_along(l))
        l[[i]] = 1
    l
}

f1 <- function(n) {
    ## bad: copy and append
    l = list()
    for (i in seq_len(n))
        l[[i]] = 1
    l
}

产生相同的结果

> identical(f0(100), f1(100))
[1] TRUE

但性能不同

> sapply(10^(1:5), function(i) system.time(f0(i))[3])
elapsed elapsed elapsed elapsed elapsed 
  0.000   0.000   0.002   0.014   0.134 
> sapply(10^(1:5), function(i) system.time(f1(i))[3])
elapsed elapsed elapsed elapsed elapsed 
  0.000   0.001   0.005   0.253  24.520

尽管这对于当前问题的规模无关紧要，但似乎应该采用更好的预分配和填充策略，以避免不得不猜测它是否相关。更好的是，*apply在这种情况下使用 orreplicate族来避免考虑它

l <- replicate(1000, matrix(rnorm(277*349), nrow=277, ncol=349), simplify=FALSE)
于 2013-06-13T00:58:01.420 回答