当数据在列表中时，我可以对代码进行矢量化吗？

Question

我正在优化我的代码，但遇到了一些问题。我知道 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 数组而不是列表。我现在要花一些时间来验证计算，以确保它们有效！

Answer 1

我的“低悬的果实”（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 或在本例中为replicate 家族，以避免不得不考虑它

l <- replicate(1000, matrix(rnorm(277*349), nrow=277, ncol=349), simplify=FALSE)